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ABSTRACT 


The  purpose  of  this  thesis  is  to  investigate  the  use  of 
Conformal  Subdomain  Basis  Functions  (CSBF)  in  the  Method  of 
Moments  (MM)  solution  of  a  thin  wire  scatterer.  The  effect  of 
using  CSBF  on  the  computed  current  and  the  scattered  field  is 
investigated  by  formulating  and  coding  the  MM  solution  for  a 
thin  wire  loop  and  comparing  the  computed  results  for  various 
loop  sizes  to  measured  data  and  two  other  KM  codes. 
Significant  reduction  in  the  number  of  segments  (and  computer 
memory  requirements)  are  found  for  loops  with  circumferences 
of  less  than  one  to  two  wavelengths  for  plane  wave  incidence. 
From  these  results,  it  is  concluded  that  the  use  of  CSBF  will 
significantly  reduce  the  number  of  segments  required  for  the 
MM  solution  of  a  spiral  antenna. 
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I.  INTRODUCTION 


Numerical  techniques  for  solving  electromagnetic  scattering  problems  using  integral 
equations  and  the  method  of  moments  (MM)  are  well  known.  The  physical  problem, 
specified  by  Maxwell’s  equations  and  boundary  conditions,  is  reduced  to  an  integro- 
differential  equation  over  finite  domains,  and  solved  using  a  procedure  referred  to  as  the 
method  of  moments  [Ref  1].  The  unknowns  (usually  currents)  are  represented 
by  a  series  of  basis  functions  with  unknown  expansion  coefficients.  The  MM  process 
generates  a  set  of  linear  equations  that  must  be  solved  simultaneously.  Until  recently, 
these  techniques  have  been  limited  to  small  (1  to  10  wavelength)  geometries  because  of 
computer  run  time  and  memory  constraints.  With  the  development  of  faster  computers 
with  more  memory,  the  MM  technique  has  increasing  application  to  larger  geometries. 
However,  computer  memory  and  run  time  can  still  be  inadequate  to  solve  many 
important  antenna  and  scattering  problems.  Numerically  efficient  solutions  require  less 
computer  memory  and/or  less  computer  run  time.  Therefore,  any  increase  in  the 
efficiency  of  a  MM  solution  is  of  great  practical  interest. 

The  usual  MM  approach  to  modeling  a  thin  wire  of  arbitrary  shape  is  to 
specify  a  series  of  points,  with  piecewise  linear  segments  between  the  points  to 
approximate  the  wire.  The  current  is  represented  by  one  or  more  basis  functions,  each 
with  constant  phase,  over  a  piecewise  linear  segment.  Typically,  the  size  of  the  segments 
is  set  by  how  accurately  the  current  or  scattered  field  needs  to  be  determined.  For 
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convergence  of  the  current,  segment  lengths  of  0.05  X  to  0.1  X  are  generally  required. 

A  second  factor  that  influences  the  segment  size  is  the  radius  of  curvature  of  the 
wire.  Tightly  curv'ed  wires  require  smaller  segments  to  reproduce  the  wire  shape 
accurately.  When  the  radius  of  curvature  is  larger  than  a  wavelength,  the  first  case  sets 
the  segment  size;  when  the  radius  of  curvature  is  much  less  than  a  wavelength,  the 
second  case  dictates  the  segment  size  (Figure  1). 

All  generally  available  MM  codes  based  on  the  method  of  subdomains  use 
the  first  approach.  A  natural  question  arises:  Does  dividing  the  wire  into  curved  segments 
that  conform  to  its  shape,  with  arclengths  restricted  only  by  the  maximum  current 
variation  rule,  yield  a  solution  that  converges  with  fewer  subsections?  If  the  answer  is 
yes,  is  the  improvement  in  convergence  worth  the  greater  complexity  and  coding  effort? 
To  resolve  this  issue,  the  following  approach  is  taken: 


Figure  1.  Left:  Piecewise  Linear  Segments,  Right:  Curved  Segments 


1. 


Formulate  the  solution  for  linearly  and  circularly  polarized  plane  wave  incidence 
for  a  geometrically  simple  shape  such  as  a  loop. 

2.  Computer  code  the  solution  in  FORTRAN. 

3.  Validate  the  solution  using  other  MM  solutions  and  measured  data. 

4.  Study  the  convergence  of  the  solution  with  respect  to  loop  parameters  and 
compare  its  performance  to  a  method  using  linear  subsections. 

5.  Study  the  effect  of  changes  in  program  structure  on  its  computational  efficiency. 
Chapter  II  discusses  the  derivation  and  the  MM  solution  of  the  electric  field  integral 

equation  for  a  thin  wire  loop.  Chapter  III  discusses  three  MM  FORTRAN  programs  used 
to  determine  the  current  on  the  loop.  The  entire  domain  solution  (due  to  R.  F. 
Harrington),  which  uses  complex  Fourier  modes  (HARLOOP)  is  considered  to  be  the 
most  accurate  and  therefore  serves  as  a  baseline  for  evaluating  the  other  solutions.  A 
second  program  that  uses  linear  segments  (LOOPSCAT)  is  compared  to  a  third  program 
that  uses  curved  segments  (CURVENEW).  Chapter  IV  discusses  the  results  obtained  by 
the  three  methods  and  presents  some  guidelines  in  choosing  an  optimum  solution  method 
for  a  given  antenna  or  scattering  problem. 


3 


n.  THE  THIN  WIRE  INTEGRAL  EQUATION 


A.  DERIV  ATION  OF  THE  THIN-WIRE  ELECTRIC  FIELD  INTEGRAL 

EQUATION 

In  this  chapter,  the  integral  equation  for  the  current  on  a  thin  wire  will  be 
developed.  Time-harmonic  field  quantities  are  assumed  throughout.  Phasor  quantities  are 
used  with  the  d"'  dependency  suppressed. 

Referring  to  the  thin  wire  geometry  of  Figure  2,  the  origin  is  point  0,  the  location 
of  a  source  point  is  given  by  the  vector  r'  and  an  observation  point  by  the  vector  r.  The 
wire  radius,  a,  is  considered  constant  over  the  length  of  the  wire.  The  vector  I  is 
everywhere  parallel  to  the  surface  of  the  wire.  Since  the  sum  of  the  tangential 
components  of  the  incident  and  scattered  electric  field  must  vanish  at  the  surface  of  a 
perfect  electric  conductor,  the  boundary  condition  is, 

ri  X  =0  (2.1) 

If  the  radius  of  the  wire  is  small  compared  to  the  wavelength  of  the  excitation,  the 
surface  current  density,  J„  can  be  considered  constant  around  the  circumference  of  the 
wire  and  directed  along  its  axis.  The  excitation  field  can  be  either  an  incident  wave  or 
an  impressed  voltage.  The  scattered  field  is  the  field  due  to  the  current  on  the  conductor 
induced  by  the  excitation  field. 

The  wave  equation  in  terms  of  the  vector  potential  A  is  given  by 
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v^a+p^a  =  -nJ^ 


(2.2) 


where  /3=2x/X.  The  solution  to  equation  (2.2)  is 


A  =  - dS'  =  ki/j,  g(ry)dS' 

AttJ  I _ 'I  J  ’ 


(2.3) 


471  |r-r- 


where  the  integration  is  over  the  primed  (source)  coordinates.  The  Green’s  function, 
g(r,r')  is  defined  as, 


g(ry)  = 


471  \r-r'\ 


(2.4) 


The  expression  for  the  scattered  electric  field  in  terms  of  A  is, 
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=  -y(jA — J—V(V-A)  . 

Ci)(ie 


(2.5) 


Applying  the  boundary  condition  of  equation  (2.1), 

=;o)A^-i-V(V-/|)  on  5  .  (2.6) 

upe 

Substitution  of  A  in  equation  (2.3)  into  equation  (2.6)  gives, 

nfjg(r,r')dS'  on  S.  (2.7) 

S' 

Call  the  second  term  on  the  right  side  of  equation  (2.7)  V,  and  assume  the  medium  to 
be  homogeneous, 

V  =  V  fv  •  p|j,  giry)dS']  =  V  fp|v  •  (J,  g(r,r'))JS'|  .  (2.8) 

.  .  s' 

The  vector  identity  for  the  divergence  of  a  scalar  u  times  a  vector  v  is. 

V  (uv)  =  Vu-v+m(V-v)  .  (2.9) 

Applying  this  identity  to  equation  (2.8)  gives 

V  =  V  pJ(Vg(r,r'))-J,d5'|  .  (2.10) 

r 

S 

It  can  be  shown  that  Vg(r,r')  =  -V’g(r,r')  [Ref.  2].  Using  this  in  equation  (2.10) 
and  applying  the  identity  of  equation  (2.9)  again  yields, 
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V  =  -vliLf(V'g(r.r'))^J,dS' 

i 

=V  L/g(r,r')(V'*J  (giry)J ,)dS' 

s'  s' 


(2.11) 


where  V  is  the  del  operator  defined  with  respect  to  the  primed  coordinates.  The  second 
integral  on  the  right  side  of  equation  (2. 1 1)  is  equal  to  zero  by  the  surface  divergence 
theorem  [Ref.  3].  Thus,  V  simplifies  to 


V  =  p/(V'-J,)V^(r,rV5'  .  (2.12) 

s 

Substitution  of  equation  (2.12)  into  equation  (2.7)  gives  an  integral  equation  for  J„. 


^'un  =  j^\^(j,8(.r,r')ds' )^g{r,r')dS'  (2.13) 

s'  s' 

which  may  be  expressed  more  compactly  as 

=  [\j^\iJ,giry)*J-{V'-J^)Vg{r,r)]dS'  .  (2.14) 

(oe  I 

s 

Equation  (2.14)  is  a  form  of  the  Electric  Field  Integral  Equation  (EFIE).  The  unknown 
quantity  to  be  solved  for  is  J,. 


B.  SOLUTION  OF  THE  EFIE  USING  MM 

The  method  of  moments  (MM)  technique  can  be  used  to  solve  for  J,  by  expanding 
it  into  a  series  of  basis  functions,  Jj, 
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(2.15) 


J,  =  EC, A 

I 

where  the  C,  are  complex  constants  to  be  determined.  Substitution  of  equation  (2. 15)  into 
2.14  gives 


B'un-'tcJ 


u>e 


dS' 


(2.16) 


To  generate  the  required  N  equations  to  solve  for  the  N  unknowns,  define  a  suitable 
weighting  function  W\.  and  take  the  inner  product  of  with  both  sides  of  equation 
2.16.  The  inner  product  is  defined  such  that  it  satisfies 


<H’,V>  =  <V,H*> 

<a/+YV,H’>  =  a  </,*♦’>  + Y  <v,H'> 
<v*,v>  >0  if  V  *  0 
<v*,v>  =0  i/  V  =  0 

[Ref.  4],  Choose  the  following  inner  product; 

<H',v>  =  jw^'vdS 
s 

where  *  signifies  the  complex  conjugate.  This  leads  to 


(2.17) 


(2.18) 


(^k-E'^dS  =  r/co^W^-52C,/[/,g(r,r')*-^(V'-J^)Vg(r,r')]<iS'dS  .  (2.19) 

1  i  <■>  i.i  J 

Interchanging  the  order  of  summation  and  integration  and  applying  the  surface  divergence 
theorem  yields  for  the  right  hand  side  of  equation  (2.19) 
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y;cJf[;(o^(W*-J.)g(r,r')--^(V'-J.)(V- •  (2-20) 

^4  J  {l)G 

1*1  c  c  ^  J 

By  making  the  following  definitions, 

K  ‘ 

.  (2.21) 

Za  =  f  f  =^(V^*J.)(V-W4)g(r,r^)  dS^dS 

s/5,L  “e 

2.20  and  2.21  can  be  written  in  matrix  form, 

[V]  =  [Z][C]  (2-22) 

where  [V].  [Z]  and  [C]  are  called  the  generalized  voltage,  impedance  and  current 

matrices,  respectively.  The  unknown  [C]  may  be  solved  by  an  appropriate  matrix 
inversion  algorithm.  Symbolically, 

[C]  =  [Zy'[V]  .  (2.23) 

The  generalized  current  matrix  elements  are  the  weighting  coefficients  in  the  summation 
of  equation  (2.15).  The  current  J,  is  computed  from  equation  (2.15)  and  the  scattered 
field  is  calculated  using  this  current  in  the  radiation  integrals.  It  should  be  noted  that  [V], 
[Z],  and  [C]  have  units  of  volts,  ohms,  and  amperes,  but  are  not  unique.  In  general,  they 
depend  on  the  choice  of  basis  and  weighting  functions.  However,  the  current  will 
converge  to  the  same  numerical  value  as  the  number  of  basis  functions  are  increased, 
provided  the  solutions  are  implemented  correctly. 
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C.  SPECIALIZATION  OF  THE  EFIE  TO  A  CIRCULAR  LOOP  USING 


CONFORMAL  SUBSECTIONS. 

The  MM  procedure  will  now  be  applied  to  a  circular  loop  in  the  X-Y  plane  as 
illustrated  in  Figure  3.  The  loop  is  an  ideal  test  geometry  to  study  the  characteristics  of 
a  MM  solution  using  conformal  subsections.  It  is  a  relatively  simple  geometry  and  other 
accurate  solution  methods  are  available  to  evaluate  the  performance  of  conformal 
subsections.  The  loop  has  a  radius  ro  and  is  divided  into  N  conformal  segments.  In  this 
case  a  conformal  segment  is  a  circular  arc.  The  arclength  of  the  i*^  segment  is, 

A',  = =  ’■M.  ■ 

The  basis  functions,  J.,  of  equation  (2.15)  are  chosen  to  be  overlapping  triangles. 


Figure  3.  Thin  Wire  Loop  Geometry 
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Triangular  basis  functions  are  chosen  because  they  are  a  more  accurate  representation  of 
the  current  than  a  pulse  basis  function  since  the  current  is  continuous  everywhere  along 
the  wire,  and  they  are  relatively  easy  to  deal  with  analytically.  Balanis  [Ref.  5] 
states  that  increasing  the  basis  function  complexity  beyond  triangles  may  not  be 
warranted  by  the  additional  improvement  in  convergence  rate.  The  triangular  basis 
functions  span  two  segments  and  overlap  as  shown  in  Figure  3.  Therefore,  the  resultant 
current  will  be  piecewise  linear.  The  weighting  (or  testing)  functions  \\\  in  equation 
(2.19)  are  chosen  such  that  (Galerkin’s  method).  Wang  [Ref.  6]  states 

that  Galerkin’s  method  provides  numerical  results  which  are  more  accurate  than  other 
testing  methods  under  similar  computational  constraints.  Substitution  of  the  above 
weighting  and  basis  functions  in  the  expression  for  in  equation  (2.21)  yields 


z^=  f  yo)p7*(/)r.(/')(/'/)--^r'.(/')r,(/)  g{ry)dS'dS 

a. 

where  t .{1')  =  -2,  T\{1)  =  . 

dl'  dl 


(2.26) 


Using  the  following  relations, 
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il'  =  4).^'  =  COS(<|) -<()') 

/  =  r(,4);  I'  =  ro(}>' 
dS  =  lizar^d^;  dS'  =  Inar^d^' 
m)  =  7’.(ro<t>)  =  r.(/) 

7'.(Z')  =  -LrV4>) 


T1  = 


g(r,r')  = 


.  -t ;  P  =  co/pi 
N  e 

7PI*I 


4n\R\ 


r-r 


equation  (2.26)  may  be  written  as 


= 


VyPn 

4n 


// 


7-*(4))r.(4)')cos(4»-4>') - ^7-/((t>')7-/((|)) 

P*'-' 


.-jO'lti 


■d^d^' 


From  Figure  4  and  the  law  of  cosines,  |  R  [  is  given  by. 


|i?|^  =  2rQ^[l  -cos(4) -(})')] *0^ 


(2.27) 


(2.28) 


(2.29) 


Figure  4.  Geometry  for  Determining  R. 


and  using  the  trigonometric  identity, 


.  In- 


sin‘ 


(1-cosa) 


(2.30) 


results  in, 


|/?l  =  /•„ 


4sm^ 


4,-4)' 


(2.31) 


By  choosing  the  test  (unprimed)  points  at  the  center  of  the  wire  and  the  source  points  on 
the  surface  of  the  wire,  the  singularities  along  the  diagonal  of  [Z]  at  <t>  =  4>'  where  r=r' 
are  avoided.  The  technique  used  to  calculate  |  R  |  is  discussed  further  in  Chapters  III 
and  IV. 

The  voltage  elements,  Vj,,  given  in  equation  (2.21)  become 


j  r*(ro4))-£‘^4>  . 


(2.32) 


The  incident  field,  E‘,  for  the  purpose  of  this  study,  will  be  a  plane  wave.  Figure  5 
shows  the  direction  of  incidence  of  the  plane  wave  in  spherical  polar  angles  6  =  6  and 
</)  =  4>  measured  from  the  Z  and  X  axes,  respectively.  E'  can  be  0  or  4>  polarized. 
Referring  to  Figure  5,  for  0  polarization,  the  component  of  E  tangential  to  the  loop,  is 


£e’e-4> 


=  EqCosQ  sin(4>-4>) 


(2.33) 


Similarly,  for  4>  polarization. 
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£^6  -4)  =  iE^cos(<I>-(t))  .  (2-34) 

The  component  of  the  phase  vector,  parallel  to  r  is 

p  -r  =  sin0  cos(<I>-(J))  .  (2.35) 

Equations  (2.30)  and  (2.32)  combine  with  2.29  to  give 

^*•5 

=  roEe' cose  J  rj(4))sm(<I>-(|))c'^'''®“®  .  (2-36) 
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A  similar  expression  for  a  <i>  directed  incident  field  is 


.  (2.37) 

♦» 

The  computer  coding  of  the  solution  for  the  thin  wire  loop  using  the  equations 
developed  above  is  described  in  the  next  chapter. 
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ffl.  COMPUTER  CODES  FOR  THE  THIN  WIRE  LOOP 


In  this  section,  the  FORTRAN  program  for  a  thin  wire  loop  using  curved  segments 
is  discussed.  The  results  are  presented  and  compared  to  similar  solutions  using  straight 
subsections  and  Fourier  modes. 

A.  DESCRIPTION  OF  THE  CODES 

Table  1  summarizes  the  organization  of  the  three  programs.  Computer  listings  are 
given  in  Appendix  B  and  equations  from  Chapter  11  will  be  referenced  with  a  "2." 
preceding  the  equation  number.  The  FORTRAN  source  code  for  the  conformal 
subsections  is  named  CURVENEW,  and  the  codes  for  the  straight  subsections  and  the 
Fourier  mode  solution  are  named  LOOPSCAT  and  HARLOOP  respectively. 

CURVENEW,  LOOPSCAT  and  HARLOOP  are  functionally  similar.  Each 
calculates  the  loop  geometry  based  on  the  segment  size,  loop  radius,  and  wire  radius,  and 
each  uses  Gaussian  quadrature  for  numerical  integration.  CURVENEW  computes  the 
impedance  matrix,  [Z],  in  subroutine  ZCURVED  from  the  loop  geometry  of  Figure  5 
using  equation  (2.28).  LOOPSCAT  uses  a  similar  formulation  applied  to  straight 
segments  in  subroutine  ZMATWW.  HARLOOP  computes  [Z]  using  the  equations 
developed  in  reference  [7],  CURVENEW  computes  the  excitation  vector,  [V],  using 
equation  (2.36)  or  (2.37)  in  subroutine  CURVEW.  LOOPSCAT  uses  a  similar 
formulation  for  straight  subsections  in  subroutine  PLANEW.  HARLOOP  computes  [V] 
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using  the  equations  in  reference  [7]  in  subroutine  PLANEW.  In  CURVENEW,  all 
integrals  are  evaluated  numerically  and  symmetry  of  the  impedance  elements  is  used  to 
fill  the  [Z]  matrix  and  reduce  the  number  of  numerical  calculations.  Two  formulations 
were  investigated  with  LOOPSCAT;  One  using  a  delta  function  approximation  to 
evaluate  one  of  the  double  integrals  in  equation  (2.21)  and  the  other  using  Gaussian 
quadrature  for  both  integrations.  CURVENEW  does  all  numerical  integrations  using 
Gaussian  quadrature.  Matrix  symmetry  is  also  used  in  LOOPSCAT  to  reduce  the  number 


TABLE  1.  FUNCTIONAL  SUMMARY  OF  PROGRAMS 


Program--  > 

Function 

.. 

CURVENEW 

(Curved 

Segments) 

LOOPSCAT 

(Straight 

Segments) 

HARLOOP 

(Fourier 

modes) 

Read  Input 
Parameters 

Lines  25-38 

Lines  15-27 

Lines  19-30 

Establish 

Loop 

Geometry 

Lines  45-71 

Lines  53-76 

Lines  44-54 

Calculate 

[Z] 

Subroutine 

ZCURVED 

Subroutine 

ZMATWW 

Subroutine 

ZMATWW 

Calculate 

[V] 

Subroutine 

CURVEW 

Subroutine 

PLANEW 

Subroutine 

PLANEW 

Solve 

System 
[C]  =  [Z]'[V] 

Subroutines 
DECOMP  and 
SOLVE 

Subroutines 
DECOMP  and 
SOLVE 

Subroutines 
DECOMP  and 
SOLVE 

Calculate  E* 

Lines  140-180 

Lines  160-197 

Lines  93-137 
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of  calculations  for  [Z].  Computation  of  [C]  is  performed  in  subroutines  DECOMP  and 
SOLVE  [Ref.  3],  which  solve  the  system  of  equations  using  Gaussian  elimination. 
Subroutines  DECOMP  and  SOLVE  are  common  to  all  three  programs.  The  subroutines 
ZCLfRVED  and  CURVEW  are  discussed  in  more  detail  in  the  next  section. 

1.  Loop  Geometry 

To  generate  the  loop  geometry  an  initial  estimate  of  the  desired  circular  arc 
length,  Al,  is  provided.  This  is  used  to  estimate  an  angular  increment,  A(t>, 

A4)  =  —  .  (3.1) 

^0 

From  this  estimate,  the  number  of  generating  points  is  calculated  by. 


N  =  Int 


+  1 


and  A6  is  recalculated  using. 


(3.2) 


A<})  =  —  (3.3) 

N 

to  ensure  that  Act)  is  such  that  N  segments  fill  exactly  2ir  radians.  The  loop  points  P,  and 
Pj  are  coincident  with  Pn  i  and  P^o  so  that  the  current  is  continuous  around  the  loop. 

2.  Subroutines  ZCITIVED  and  ZMATWW 

Subroutines  ZCURVED  and  ZMATWW  take  advantage  of  the  symmetry  that 
exists  on  the  loop.  For  all  basis  functions,  the  self  impiedance  terms  are  equal.  The 
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remainder  of  the  matrix  is  filled  with  impedance  values  that  repeat  in  a  cyclic  manner. 
Mathematically, 


Zn 

~  ^33  "• 

II 

N 

r  _ 

2  ,  = 

Z,  = 

=  2 

'12 

^23 

■^34  "• 

S 

r  _ 

2  = 

2,,  = 

=  2 

'13 

■^24 

■^35  "• 

^N-2  N 

(3.4) 


2  =2 

The  elements  along  any  diagonal  of  [Z]  are  equal  and  the  lower  off-diagonal  elements 
are  the  mirror  image  of  the  upper  diagonals.  Thus  [Z]  is  a  symmetrical  Toeplitz  matrix. 
Therefore,  compulation  of  the  first  row  of  [Z]  provides  enough  information  to  fill  the 
entire  matrix. 

Because  of  the  Green’s  function  in  the  integrand  for  the  impedance  elements,  the 
numerical  treatment  of  the  self  term  is  very  important.  To  optimize  the  convergence  rate 
and  accuracy  of  CURVENEW  and  LOOPSCAT,  several  different  approaches  are  used 
to  evaluate  |  R  j  near  the  singularity  point  where  r  =  r'.  In  the  first  method,  the 
obser\ation  point  is  chosen  along  the  axis  of  the  wire  and  the  source  point  along  the 
surface  for  all  i,k,  giving  |  R  1  =a  at  <i)=<i>'  (equation  (2.31)).  For  the  second  method, 
both  the  obser\ation  point  and  the  source  point  are  chosen  along  the  axis  of  the  wire 
except  on  the  segment  i=k  where  the  value  of  at  the  midpoint  is  chosen  on  the  axis 
of  the  wire,  with  r'  on  the  surface  of  the  wire.  Finally,  both  the  observation  point  and 
the  source  point  are  chosen  along  the  axis  of  the  wire,  except  on  the  segment  i  =  k,  where 
r  is  chosen  along  the  axis,  and  r'  is  chosen  on  the  surface.  Choosing  the  source  point  and 
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observation  point  as  in  the  first  case  gives  the  most  accurate  results,  but  only  slightly 
more  accurate  than  the  third  case.  Case  two  is  accurate  for  small  segment  sizes  but  is 
inaccurate  for  larger  segment  sizes.  Case  three  was  selected  for  both  CURVENEW  and 
LOOPSCAT  because  it  is  only  slightly  less  accurate  than  case  one,  and  required  fewer 
lines  of  code. 

ZCURVED  calculates  the  impedance  elements  of  the  first  row  of  [Z]  by  breaking 
the  integral  in  equation  (2.28)  into  four  parts.  For  example,  in  the  first  row  of  [Z],  Z,,, 
is  calculated  by  summing  contributions  from  the  following  four  regions  of  integration 
(Figure  6): 

1 .  A  double  integration  along  the  positive  slope  of  the  T,  over  segment  1  and  the 
positive  slope  of  T,  over  segment  i. 

2.  A  double  integration  along  the  negative  slope  of  the  T,  over  segment  2,  and  the 
positive  slope  of  T,  over  segment  i. 

3.  A  double  integration  along  the  positive  slope  of  the  Tj  over  segment  1,  and  the 
negative  slope  of  T,  over  segment  i  + 1. 

4.  A  double  integration  along  the  negative  slope  of  the  T,  over  segment  2.  and  the 
negative  slope  of  T,  over  segment  i  +  1. 

The  integrations  are  computed  in  a  similar  manner  for  the  derivatives  of  the  basis 
functions,  T,'  and  T,',  over  the  same  subsections.  A  similar  procedure  is  used  for  the 
straight  subsections  in  subroutine  ZMATWW'  to  calculate  the  impedance  elements. 

The  numerical  integrations  are  performed  using  Gaussian  quadrature,  with 
the  number  of  points  per  subsection  specified  as  an  input  parameter  to  program 
GAUSWGT,  which  computes  the  Legendre  polynomial  coefficients  for  a  sjjecified 
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number  of  integration  points  and  writes  them  to  file  OUTGLEG.  Gaussian  quadrature 
was  chosen  because  it  requires  fewer  function  evaluations  than  other  methods  for  a  given 
accuracy  and  does  not  require  equal  interval  samples  [Ref.  8].  CURVENEW, 
LOOPSCAT  and  HARLOOP  read  the  coefficients  from  file  OUTGLEG.  The  number 
of  integration  points  per  wavelength  was  varied  to  optimize  convergence  of  LOOPSCAT 
and  CURVENEW  and  is  discussed  in  more  detail  in  Chapter  IV.  The  excitation  vector 
[V]  is  calculated  from  equation  (2.36)  or  (2.37)  in  subroutine  CURVEW  of 
CURVENEW. 

3.  Execution  Time 

Analysis  of  the  nested  DO  loop  structure  of  subroutine  ZCURVED  of 
program  CURVENEW  indicates  that  the  total  execution  time  of  ZCURVED  can  be 
represented  by. 

T.  =  a  N 

aC  ^  I  ^ 


Figure  6.  Loop  Geometry  for  Impedance  Integrations 
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where  is  the  number  of  curved  segments  (from  equation  (3.2))  is  the  number  of 
Gaussian  integration  constants  per  curved  segment  and  is  a  constant.  Execution  time 
of  the  subroutines  DECOMP  and  SOLVE,  common  to  both  CURVENEW  and 
LOOPSCAT,  can  be  represented  similarly  by  [Ref.  9], 

Ti  = 

where  N  is  the  number  of  segments.  The  excitation  subroutine  CURVEW  execution  time 
and  field  integrations  are  given  by, 

T,,  =  (3.7) 

where  again  y  and  are  constants.  Assuming  that  the  execution  time  of  the  rest  of  the 
program  is  negligible,  the  total  execution  time  of  CURVENEW  is 


(3.8) 


A  similar  expression  for  LOOPSCAT  uses  the  subscript  1, 


•  (3.9) 

Run  times  were  recorded  for  various  values  of  N^,  N,  and  using  an  IBM  PC/AT  with 
a  math  coprocessor  and  the  coefficients  for  CURVENEW  are  found  to  be  7=0.000156, 
Qfj =0.0230,  and  ^,=0.0222.  The  coefficients  for  LOOPSCAT  are  ai=0.0132  and 
^,=0.0252.  For  the  moment,  assume  that  the  number  of  Gaussian  integration  points  per 
wavelength,  Ng,  is  held  constant  for  both  CURVENEW  and  LOOPSCAT.  The  number 
of  integration  points  on  a  segment  is. 
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and  the  number  of  segments  is, 


(3.11) 


Similar  expressions  may  be  wntitr.  ^or  straight  subsections.  Combining  equations  (3.8), 
3.10,  and  3.11  gives 


g 


(3.12) 


The  ratio  of  to  T,  is  given  by. 


An^rlNlaJ  +2TirQlic^g 

47tVojVga,/iV^  +  YiV/  +  27iro(r^iV^ 


TJT,  = 


(3.13) 


Equation  (3.13)  will  be  used  in  Chapter  IV  to  compare  the  execution  times  of 


CURVENEW  and  LOOPSCAT. 


rV.  CALCULATED  DATA  FOR  THE  LOOP 


The  convergence  of  the  MM  solutions  for  both  the  current  and  electric  field  for 
circular  loops  of  various  dimensions  are  presented  for  both  linear  and  circular 
polarizations.  The  convergence  is  shown  to  depend  on  the  segment  size  and  number  of 
integration  points,  as  well  as  excitation  conditions  (incidence  direction  and  polarization). 
Representative  plots  are  presented  within  the  chapter,  and  additional  plots  are  given  in 
Appendix  B. 

A.  CONVERGENCE  OF  HARLOOP 

The  Fourier  mode  solution,  HARLOOP,  was  tested  for  convergence  with  respect 
to  incidence  angle,  number  of  modes,  and  number  of  integration  constants  to  establish 
a  baseline  for  comparison  to  CURVENEW  and  LOOPSCAT.  HARLOOP  was  chosen  as 
a  baseline  because  the  sinusoidal  basis  functions  match  the  physical  behavior  of  the 
current  on  the  loop,  and  thus  the  current  series  converges  rapidly.  This  is  illustrated  in 
Figure  7  for  a  0.5  X  radius  loop  with  a  plane  wave  incident  at  an  angle  of  40  degrees. 

Oscillation  of  the  current  as  a  function  of  4>  becomes  more  rapid  as  0  is  increased, 
because  the  phase  of  the  incident  field  over  the  loop  varies  as  sin  0  (equation  (2.35)). 
For  a  0  polarized  linear  wave  incident  in  the  </>=0  plane,  the  current  is  always  zero  at 
<t>=0  and  180  degrees,  where  E'  is  cross-polarized  with  respect  to  the  axis  of  the  wire. 
For  6  polarized  incident  waves,  maxima  occur  at  <i>=90  and  270  degrees,  where  E'  is 
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X  1 


Currar^t  Mogrtituda  vs  Uoop  ongia 


Figure  7.  Convergence  of  the  Current  in  the  Complex  Exponential  Solution 

parallel  to  the  axis  of  the  wire.  The  overall  current  amplitude  decreases  for  0 
approaching  90  degrees,  as  expected,  since  the  loop’s  projected  area  is  small  as  viewed 
by  the  incident  wave.  For  polarized  incident  waves,  the  minima  occur  at  <^>=90  and 
270  degrees  and  maxima  at  <^>=0  and  180  degrees.  The  currents  do  not  vanish  for  0 
approaching  90  degrees  because  the  loop  is  parallel  to  the  </>  polarized  incident  field. 
Circularly  polarized  incident  waves  give  a  constant  magnitude  current  at  normal 
incidence,  and  oscillations  increase  with  0.  For  0  =  90  degrees,  the  circular  and  <t> 
polarization  responses  are  identical. 

HARLOOP  is  also  found  to  be  in  agreement  with  measurements  taken  on  the  echo 
area  of  wire  loops  at  normal  incidence  [Ref.  10].  The  plot  of  Figure  8  gives  the 
echo  area  (a/X^)  versus  Tq  for  varying  wire  radius  using  HARLOOP.  Measured  data  is 
indicated  by  the  ’  +  ’  sign. 
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B.  CONVERGENCE  OF  THE  CURRENT  EXPANSION 


Having  established  HARLOOP  as  a  baseline  for  comparison,  convergence  of  the 
curved  subsection  program  CURVENEW  and  the  linear  subsection  program  LOOPSCAT 
was  evaluated.  The  plots  of  Figures  9  through  14  give  the  current  on  the  loop  as  a 
function  of  loop  angle,  <t>,  for  varying  0,  loop  radius,  and  incident  wave  polarization. 
The  number  of  integration  points  per  wavelength,  N,,  is  held  constant  at  320  in 
LOOPSCAT  and  CURVENEW.  This  number  was  chosen  empirically  to  give  a 
converged  current  within  five  to  ten  percent  RMS.  The  RMS  error  is  defined  relative  to 
the  Fourier  mode  solution.  The  HARLOOP  current  is  plotted  with  the  solid  line,  those 
of  CURVENEW  are  plotted  with  the  "  +  "  sign,  and  those  of  LOOPSCAT  with  the  "o". 
The  wire  radius  is  0.001  X  for  these  calculations. 


Bockacott*''  Ec^^o  v»  Looo  Podlw* 


Figure  8.  Backscatter  Echo  Area  for  a  Loop  with  varying  Radius  at  Normal  Incidence 
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Figure  9.  Magnitude  of  the  Current  on  a  0. 1  X  Radius  Loop,  Normal 
Incidence,  Circular  Polarization  (+  =curved;  o  =linear) 
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Figure  10.  Phase  of  the  Current  on  a  0.1  X  Radius  Loop,  Normal 
Incidence,  Circular  Polarization  (+  =curved;  o  =  linear) 


27 


xio- 


CufTtnt  Mcgnitud*  V*  Loop  onglo 


Loop  Anglo,  dogroot 

Figure  11.  Magnitude  of  the  Current  on  a  0.1  X  Radius  Loop,  Incidence 
Angle=40  deg,  Circular  Polarization  (+  =curved;  o  =linear) 
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Figure  12.  Phase  of  the  Current  on  a  0.1  X  Radius  Loop,  Incidence 
Angle =40  deg..  Circular  Polarization  (+  =  curved;  o  =  linear) 
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Figure  13.  Magnitude  of  the  Current  on  a  O.I  X  Radius  Loop,  Incidence 
Angle=40  deg..  Theta  Polarization  (+  =curved;  o  =linear) 
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Figure  14.  Phase  of  the  Current  on  a  0.5  X  Radius  Loop,  Incidence 
Angle=40  deg..  Theta  Polarization  (+  =curved;  o  =linear) 
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Representative  plots  of  the  normalized  root  mean  squared  error  are  given  in  Figures  15 
through  22.  For  set  values  of  N  and  Ng,  CURVENEW  converges  faster  than 
LOOPSCAT  in  most  cases.  The  error  difference  is  most  pronounced  for  loop 
circumferences  on  the  order  of  a  wavelength  or  less.  From  the  plots,  for  ro=0.1  X, 
CURVENEW  converges  to  less  than  10  percent  error  for  segment  sizes  ranging  from 
0.02  to  0.2  wavelengths  (N,=32  to  N<.=3).  LOOPSCAT  converges  to  within  10  percent 
error  for  segment  sizes  less  than  approximately  0.06  X  (N,  >  10)  but  gives  errors  of  30 
to  40  percent  for  a  segment  size  of  0.2  wavelengths.  There  is  no  improvement  using 
curved  subsections  on  larger  loops  for  off-axis  incidence  waves,  but  the  cur\'ed 
subsections  give  small  improvements  for  large  loops  at  normal  incidence.  This  is 
expected  in  view  of  the  behavior  of  the  current  on  the  loop.  For  linear  polarization  the 
current  abruptly  flips  polarity  from  one  side  of  the  loop  to  the  other.  For  circular 
polarization,  the  amplitude  is  constant,  but  the  phase  is  linear.  Both  of  these  conditions 
can  be  represented  accurately  by  a  few  triangles  if  the  impedance  and  excitation  integrals 
are  evaluated  precisely  on  the  loop  contour  (see  Figure  23). 

Plots  of  the  magnitude  of  the  backscattered  E'  versus  incidence  angle  for  varying 
segment  size,  Ng,  and  ro  are  given  in  Figures  24  and  25.  As  with  the  currents,  the 
backscattered  field  converges  more  rapidly  for  CURVENEW  than  LOOPSCAT  in  most 
cases,  with  the  greatest  difference  for  small  radius  loops  and  angles  near  the  maximum 
values  of  |  E’  ]  . 
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Figure  15.  Error  in  the  Current  Magnitude  for  a  0.1  X  Radius  Loop, 
Normal  Incidence,  Circular  Polarization  (+  =curved;  o  =linear) 


Figure  16.  Error  in  the  Current  Magnitude  for  a  0. 1  X  Radius  Loop, 
Incidence  Angle=40  deg..  Circular  Polarization  (+  =curved;  o  =linear) 
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Figure  17.  Error  in  the  Current  Magnitude  for  a  0.5  X  Radius  Loop, 
Normal  Incidence,  Circular  Polarization  (+  =curved;  o  =Iinear) 
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Figure  18.  Error  in  the  Current  Magnitude  for  a  0.5  X  Radius  Loop, 
Incidence  Angle =40  deg..  Circular  Polarization  (+  =curved;  o  =linear) 
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Figure  19.  Error  in  the  Current  Magnitude  for  a  0.1  X  Radius  Loop, 
Normal  Incidence,  Theta  Polarization  (+  =curved;  o  =linear) 
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Figure  20.  Error  in  the  Current  Magnitude  for  a  0.1  X  Radius  Loop, 
Angle  of  Incidence=40  deg..  Theta  Polarization  (+  =curved;  o  =linear) 
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Figure  21.  Error  in  the  Current  Magnitude  for  a  0.5  X  Radius  Loop, 
Normal  Incidence,  Theta  Polarization  (+  =curved;  o  =linear) 
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Figure  22.  Error  in  the  Current  Magnitude  for  a  0.5  X  Radius  Loop, 
Incidence  Angle =40  deg..  Theta  Polarization  (+  =cur\'ed;  o  =  linear) 
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Figure  23.  A  Representation  of  a  Sinusoidal  Current  with  Two  Basis  Functions  (top)  and 
a  Constant  Current  as  a  Superposition  of  Triangles  (bottom) 
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Figure  24.  Backscattered  Electric  Field  Intensity  for  varying  Angles  of 
Incidence,  0. 1  X  Radius  Loop,  Theta  Polarization  (+  =curved;  o  =  linear) 
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Figure  25.  Backscattered  Electric  Field  Intensity  for  varying  Angles  of 
Incidence,  0.5  X  Radius  Loop,  Theta  Polarization  (+  =curved;  o  =linear) 
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C .  COMPARISON  OF  EXECUTION  TIME  ANT)  MEMORY  REQUIREMENTS 
The  average  run  time  for  a  given  loop  radius  and  convergence  error  is  greater  for 
CURVENEW  than  for  LOOPSCAT  due  to  the  Ng^  dependence  in  equation  (3.13)  and 
relative  magnitudes  of  the  coefficients  a  and  7.  The  plot  of  equation  (3.13)  in  Figure  26 
illustrates  the  ratio  of  run  times  of  CURVENEW  and  HARLOOP  versus  N,  for  Nc=4, 


Retie  of  R^n  Tirnes  vs.  Nl 


Figure  26.  Ratio  of  Execution  Times  for  CURVENEW  and  LOOPSCAT  for  Fixed 
Curved  Segment  Lengths  and  with  Varying  Linear  Segment  Lengths 

8,  32  and  Ng  =  320  for  a  0.1  X  loop.  The  run  time  of  CURVENEW  is  less  than 

HARLOOP  for  T/T,  <  1.  From  the  plot,  the  "break  even"  points  are  approximately 

N,=  115,  90,  52  for  N,=4,  8,  32.  For  N,  less  than  these  values,  the  integration 

subroutine  ZCURVED  is  the  determining  factor  in  the  run  time;  for  N,  greater  than  these 

values,  Gausssian  elimination  subroutines  DECOMP  and  SOLVE  are  the  determining 

factors.  The  increased  number  of  integrations  per  segment  completely  offset  the  time 
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savings  of  a  reduced  [Z]  matrix  in  CURVENEW.  As  mentioned  in  Chapter  III,  a  delta 
function  approximation  for  the  outer  integration  of  the  impedance  integral  of  equation 
(2.28)  was  investigated.  This  reduces  the  exponent  of  to  one  in  equation  (3.13),  but 
many  more  segments  are  required  for  a  given  accuracy.  A  more  efficient  integration 
scheme  using  a  large  number  of  integration  points  in  the  vicinity  of  </>=</>’  and  fewer 
integrations  elsewhere  may  reduce  the  execution  time. 

The  savings  in  computer  memory  is  significant  for  CURVENEW,  since  the  number 
of  matrix  elements  is  on  the  order  of  For  a  given  accuracy  for  a  0. 1  X  loop,  the  ratio 
of  the  number  of  elements  required  for  curved  and  linear  subsections,  (N,/Nl)^  is  on  the 
order  of  0.1  to  0.2.  This  is  a  reduction  of  80  to  90  percent.  Although  the  memory 
requirements  for  the  small  loops  considered  here  are  not  prohibitive  for  piecewise  linear 
segments,  greater  memory  savings  will  be  realized  for  larger  geometries  comprised  of 
many  small  features. 
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V.  CONCLUSIONS 


The  use  of  conformal  subdomain  basis  functions  (curved  subsections)  to  represent 
the  current  on  a  thin  curved  wire  was  investigated  by  solving  the  thin  wire  electric  field 
integral  equation  using  the  method  of  moments.  A  solution  using  triangular  basis 
functions  was  computer  coded  in  FORTRAN  and  validated  by  comparing  it  to  measured 
data  and  the  results  of  two  other  method  of  moments  solutions  (LOOPSCAT  and 
HARLOOP).  The  effect  of  varying  loop  radius,  segment  size,  number  of  integration 
points  and  incident  wave  parameters  on  the  accuracy  and  rate  of  convergence  of  the 
current  expansion  and  backscattered  field  was  investigated. 

For  small  loops  with  circumferences  on  the  order  of  a  wavelength,  the  number  of 
segments  required  to  converge  to  a  given  accuracy  with  the  cursed  segments  was  as 
small  as  20  percent  of  the  number  of  linear  segments  required  to  converge  to  the  same 
accuracy  (see  Table  2).  From  computed  data,  it  was  determined  that  the  greatest 
reduction  in  the  number  of  unknowns  for  curved  subsections  occurs  for  geometries  where 
the  current  amplitude  variations  over  the  surface  are  small  and  the  phase  variations  are 
small  or  linear.  As  mentioned  in  Chapter  I,  the  general  rule  of  thumb  for  the  length  of 
one  segment  is  0.05  X  to  0.1  X,  which  corresponds  to  a  phase  variation  of  20  to  40 
degrees.  This  restriction  in  phase  variation  is  the  driving  factor  when  choosing  the 
segment  size  for  geometrical  features  having  radii  of  curvature  of  approximately  X/2  or 
larger.  A  piecewise  linear  segment  is  small  enough  to  represent  the  geometry  accurately 
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in  this  case.  The  phase  restriction  applies  to  curved  segments  as  well,  but  the  curved 
segments  conform  exactly  to  the  wire,  and  hence  for  small  loops  there  is  no  sacrifice  in 
geometrical  accuracy  by  choosing  segment  sizes  of  approximately  0.05  X  or  20  electrical 
degrees.  As  the  loop  becomes  larger,  or  the  wavelength  becomes  smaller,  the  curved  and 
straight  subsection  solutions  become  equivalent.  The  greatest  advantage  in  using  curved 
subsections  to  reduce  the  number  of  segments  is  for  electrically  small  structures  where 
small  linear  segments  are  required  simply  to  reproduce  the  wire  shape. 

Although  the  number  of  segments  was  greatly  reduced  using  conformal  subsections, 
the  execution  time  was  increased  due  to  the  increased  number  of  integration  points  per 
segment  required  for  acceptable  accuracy.  To  reduce  the  integration  time,  it  is  suggested 
that  the  number  of  integration  points  per  wavelength  be  varied  from  a  large  number  when 
evaluating  the  self  term,  to  fewer  points  away  from  the  self  term.  For  certain  geometries, 
symmetry  could  also  be  used  to  reduce  the  integration  time. 

To  avoid  singularities,  the  MM  testing  procedure  was  performed  along  the  axis  of 
the  wire  and  the  current  constrained  to  the  surface  of  the  wire.  A  delta  function 
approximation  in  the  impedance  integrations  was  found  to  reduce  the  time  required  to 
compute  the  impedance  matrix,  but  as  expected,  required  more  segments  for  a  given 
convergence  accuracy. 

A  disadvantage  in  the  formulation  of  CURVENEW  is  that  an  analytic  expression 
for  the  curve  is  needed  to  perform  the  integrations  for  the  impedance  matrix  and  the 
excitation  vector.  The  expressions  will  change  each  time  a  new  curve  is  analyzed,  and 
consequently,  considerable  effort  W'ill  be  required  to  modify  the  code  every  time  a  change 
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is  made  in  the  geometry.  Programs  that  use  linear  segments  are  more  flexible  because 
they  require  only  the  coordinates  of  the  points  along  the  wire  axis  to  generate  the 
integration  points. 

The  next  logical  step  in  testing  the  effectiveness  of  conformal  subdomain  basis 
functions  is  to  formulate  the  solution  for  an  equiangular  spiral  wire.  The  equiangular 
spiral  is  used  in  broadband  antennas  and  has  a  simple  mathematical  form.  For 
geometrical  accuracy,  the  segment  size  in  the  piecewise  linear  formulation  will  be  much 
smaller  than  20  electrical  degrees  near  the  center  of  the  spiral.  Equal  length  conformal 


TABLE  2.  COMPARISON  OF  CURVENEW  AND  LOOPSCAT 


Curved 

Subsections 

Linear 

Subsections 

Loop  Radius 

O.IX 

0.5X 

O.IX 

0.5X 

Number  of  Segments  for  10% 
RMS  Error.  Normal  Incidence 

3 

5 

21 

5 

Number  of  Segments  for  10% 
RMS  Error,  Off  Axis  Incidence 

3 

26 

10 

26 

Execution  Time*,  10%  RMS 
Error,  Normal  Incidence 

314  s 

4669  s 

32  s 

m 

Execution  Time*,  10%  RMS 
Error,  Off  Axis  Incidence 

314  s 

918  s 

59  s 

540  s 

[Z]  Matrix  Size,  10  %  RMS 
Error,  Normal  Incidence 

9 

25 

441 

25 

[Z]  Matrix  Size.  10  %  RMS 
Error,  Off  Axis  Incidence 

9 

676 

100 

676 

*  Execution  time  measured  with  an  IBM  PC/AT. 
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segments  may  be  used  along  the  spiral  arms,  and  it  is  anticipated  that  the  number  of 
required  segments  will  be  substantially  reduced. 
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APPENDIX  A 


COMPUTER  CODES 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 


C  MAIN  PROGRAM  *CURVNEW.FOR* 

C  ****  MORE  NUMERICALLY  EFFICIENT  THAN  CURVSUB.F  •*** 

C  PLANE  WAVE  SCATTERING  FROM  A  CIRCULAR  LOOP  IN  THE  Z  PLANE. 

C  METHOD  OF  MOMENTS  WITH  CURVED  SUBSECTIONS 
C  ***♦  RADAR  CROSS  SECTION  CALCULATION 
C 

COMPLEX  Z(25000),B(500),C(500),BP(500),CP(500),EX,EC 

COMPLEX  ET.EP,UC,RW(500).U0 

DIMENSION  ECP(400),IPS(500).ANG(400).EDP(400) 

DIMENSION  XG(128).AGfl28).PH(500).DEL(500).PA(500) 

DIMENSION  ECV(500),EXV(500),PHC(500).PHX(500) 

DATA  Pl/3. 14159265/ 

DATA  IPRlNT/l/.ITEST/i; 

RAD=P1/180. 

ECX  =  0. 

BK  =  2.*PI 
ETA  =  377. 

U0  =  (0.,0.) 

UC  =  (0.,-l.)*ETA*BK/4./Pl 
C0NST  =  16.*P1**3 
PHIRD  =  0. 

C 

C  READ  INPUT  AND  PROGRAM  CONTROL  PARAMETERS 
C 

OPEN(l.FILE='PARAMLST.DAT’) 

READd, *)  ANGLE. DT.START.STOP.AW.RO.SEG.HARR.GCONST.XLOC.XIPOL 
LOC  =  INT(XLOC) 

IPOL=INT(XIPOL) 

CLOSE(l) 

C 

C  READ  GAUSSIAN  CONSTANTS 
C 

OPEN(2.FILE=  ’OUTGLEG’) 

READ(2,*)  NG 
DO  I  1=1. NG 
READ(2.*)  XG(1).AG(I) 

1  CONTINUE 
C 

C  OUTCURV  IS  THE  FILE  THAT  THE  SCATTERED  FIELD  DATA  IS  WRITTEN  TO 
C  ICURVOUT  IS  THE  FILE  THAT  THE  LOOP  CURRENT  DATA  IS  WRITTEN  TO 
C 
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OPEN(8,FILE  = ’OURCURV.DAT') 
0PEN(7,FILE=  ’ICURVOUT.DAT’) 


/ 
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44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 
7*' 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 


C 

C  GENERATE  THE  LOOP  POINTS 

C 

C 

C  CHOOSE  THE  NUMBER  OF  POINTS  BASED  ON  THE  VALUE  OF  SEG 
C 

C  AW  =  AW*R0 
DPH1  =  SEG/R0 
NP=INT(2.*PI/DPH1)  +  1 
DPHI  =  360.  /FLOAT(NP) 

PH(1)  =  0. 

DO  10  1  =  2,NP+1 
PHfl)  =  FLO  AT(I- 1  )*DPHI*RAD 
DEL(1-1)  =  (PH(1)-PH(I-1)) 

PA(I-l)  =  (PH(I)  +  PH(I-l))/2. 

10  CONTINUE 
NP=NP+2 
C 

C  OVERLAP  THE  ENDS  SO  TH.AT  CURRENT  WILL  BE  CONTINUOUS  ON  THE  LOOP 
C 

PH1NP)  =  BK*PH(2) 

DEL(NP-1)  =  DEL(1) 

P.A(NP-1)  =  BK  +  PA(1) 

MT=NP-2 
DO  52  1=1,NP 
XHB  =  R0*COS<PH(I)) 

YHB  =  R0*S1N(PH(I)) 

52  CONTINUE 

WRITElO,*)  'GEOMETRY  DEFINED' 

IF(ITEST.EQ.O)  GO  TO  98 
C 

C  DEFINE  DIMENSIONS  OF  THE  IMPEDANCE  MATRIX  BLOCKS 
C 

WRITE(6,*)  'NP,MT=',NP,MT 
C 

C  COMPUTE  IMPEDANCE  MATRIX  ELEMENTS 
C 

CALL  ZCURVED(NP.RO,PH.DEL,PA.NG.XG,AG,AW,Z,ZMAX) 

DO  11  1  =  1, MT 
CZ  =  CABS(Z(I)) 

AZ  =  ATAN2(  AIM  AG(Z(I)).REAL{Z(1))  +  1  .E-8)/RAD 
1 1  CONTINUE 

>\TIITE(6.*)  'WIRE  IMPEDANCE  COMPUTED' 

C 

C  PERFORM  LU  DECOMPOSITION 
C 

CALL  DECOMP(MT.IPS.Z) 

WRITE(6.*)  'Z  DECOMPOSED' 
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93 

94 

95 

96 

97 

98 

99 

100 
101 
102 

103 

104 

105 

106 

107 

108 
109 

no 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 


C 

C  BEGIN  FIELD  CALCULATIONS 
C 

98  PHRO=PHIRD*RAD 
rr=INT((STOP-START)/DT)+ 1 
WRITE(7,*)  IT.MT.O.O 

DO  5001=  LIT 

THETA = FLO  AT(I- 1  )*DT + START 

THR=THETA*RAD 

PHR=PHR0 

IFCTHETA.LT.  180.)  GO  TO  99 
THR  =  (360.  -THET  A)*RAD 
PHR  =  PHR0  +  PI 

99  CONTINUE 
ET=U0 

EP  =  U0 
C 

C  COMPUTE  THE  EXCITATION  VECTOR 
C 

CALL  CURVEW{NP.RO,PH,DEL.PA.NG,XG,AG,THR.PHR.RW) 
IF  (LOC  .EQ.  Oj  THEN 
C 

C  CIRCULAR  POLARIZATION  IF  LOC=l  ELSE  LINEAR 
C 

IF(IPOL.EQ.l)THEN 

C 

C  THETA  POLARIZED  INCIDENT  WAVE  (IPOL=  1) 

C 

DO  101  L=1.MT 
101  B(L)  =  RW(L) 

ELSE 

C 

C  PHI  POLARIZED  INCIDENT  WAVE  (IPOL  =  2) 

C 

ENDIF 

1F(ITEST,EQ.0)  GO  TO  9998 
C 

C  PERFORM  GAUSSIAN  ELIMINATION  TO  DETERMINE  [C] 

C 

CALL  SOLVE(MT.lPS.Z,B.C) 

DO  210  L=1,MT 

WRITE(7.*)  L.L*DPHLCABS(C(L)),ATAN2(REAL(C(L)), 

*  AIMAG(C(L)))/RAD 
ET  =  ET  +  RW(L)*C(L) 

210  EP  =  EP  +  RW(L  +  MT)*C(L) 

ELSE 

C 

C  THETA  POLARIZED  INCIDENT  WAVE 
C 

DO  221  L=1.MT 
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143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 


221  B(L)  =  RW(L) 

C 

C  PHI  POLARIZED  INCIDENT  WAVE 
C  PHASE  SHIFT  FOR  CP  IS  PI/2. 

C 

DO  222  L=1,MT 

222  BP(L) = RW(L + MT)*CEXP(CMPLX(0.  .PI/2.)) 

CALL  SOLVE(MT,IPS,Z,B,C) 

CALL  SOLVE(MT,IPS.Z,BP,CP) 

DO  220  L=1.MT 

WRITE(7.*)  L.L*DPHI,CABS(C(L)  +  CP{L)),ATAN2(REAL(C(L)  +  CP(L)), 
*  AIMAG(C(L)  +  CP(L)))/RAD 

ET = ET + RW(L)*C(L)  +  RW(L)*CP(L) 

220  EP  =  EP  +  R  W(L  +  MT)*C(L)  +  RW(L  +  MT)*CP(L) 

ENDIF 
EC  =  ET*UC 
EX  =  EP*UC 
ANG(I)  =  THETA 
ECV(I)  =  CABS(EC) 

EXV(I)  =  CABS(EX) 

ECR  =  REAL(EC1 
ECI  =  AIMAG(EC) 

EXR  =  REAL(EX) 

EXI  =  AIMAG(EX) 

PHC(I)  =  ATAN2(ECLECR+1.E-20)'RAD 
PHXd)  =  ATAN2(EXI.EXR  +  1  .E-20)/RAD 
ECX  =  AMAX1(ECX.ECV(I),EXV(I)) 

500  CONTINUE 

\MIITE(6.*)  ’EMAX  =  ’.ECX 
DO  600  I=1.IT 

ECV(I)  =  AMAX1(ECV(I).1.E-10) 

EXV(I)  =  AMAXI(EXV{I),1.E-10) 

ECP(I)  =  (ECV(I)/ECX)**2 
EDP(I)  =  {EXV(I)/ECX)**2 
ECP(I)  =  AMAX1(ECP(I).. 00001) 

EDP(  I )  =  AM  AX  1  (EDPd ),  .0000 1 ) 

ECP(I)=  10.*ALOG10(ECPd)) 

EDP(I)=10.*ALOG10(EDP(I)) 

600  CONTINUE 

SIGMA  =  (ECX**2)*CABS(UC)/(2.*BK) 

SIGDB=  10.*ALOG10(S1GMA) 

WRITE(6,*)  ’BACKSCATTER  CROSS-SECTION,  IN  DB  =  ’, SIGMA. SIGDB 
208  FORMAT(/.5X.’SIGMAAVAVL  SQ  =  ’,E15.4. 

♦/,5X,'  INDB=\F8.4) 

DO  9000  L=1.IT 
WR1TE(8.*)  ANG(L).ECV(L) 

9000  CONTINUE 
9998  STOP 
END 

SUBROUTINE  SOLVE(N,IPS.UL,B.X) 


47 


193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 
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221 
222 

223 

224 

225 

226 

227 

228 
'|'^g 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 


C 

C 

C  SUBROUTINE  TO  SOLVE  SYSTEM  OF  EQUATIONS  WITH  COMPLEX 
COEFFICIENTS. 

C  CALL  ’DECOMP’  FIRST.  (FROM  MAUTZ  AND  HARRINGTON) 

C 

COMPLEX  UL(50000).B(200),X(200),SUM 
DIMENSION  IPS(500) 

NP1=N+1 
IP  =  IPS(1) 

X(1)  =  B(IP) 

DO  2  1  =  2, N 
IP  =  IPS(I) 

IPB  =  IP 
IM1  =  I-1 
SUM=(0.,0.) 

DO  1  J  =  1,IM1 

SUM  =  SUM  +  UL(IP)*X(J) 

1  IP  =  IP  +  N 

2  X{I)  =  B(IPB)-SUM 
K2  =  N*(N-1) 

IP  =  IPS(N)  +  K2 
X(N)  =  X(N)/UL(IP) 

DO  4  IBACK  =  2,N 

I  =  NP1-IBACK 

K2  =  K2-N 

IPI  =  IPS(I)  +  K2 

IP1=I+1 

SUM=(0.,0,) 

IP  =  IPI 

DO  3  J  =  IP1.N 
IP  =  IP  +  N 

3  SUM  =  SUM  +  UL(IP)*X(J) 

4  X(I)  =  (X(I)-SUM)/UL(IPI) 

RETURN 

END 

SUBROUTINE  DECOMP(N.IPS.UL) 

C 

C  SUBROUTINE  TO  DECOMPOSE  SYSTEM  OF  EQUATIONS. 

C  FROM  MAUTZ  AND  HARRINGTON. 

C 

COMPLEX  UL(50000).PIVOT,EM 
DIMENSION  SCL(200),IPS(200) 

DO  5  I=I,N 
IPS(I)  =  I 
RN  =  0. 

J1  =  I 

DO  2  J=I.N 

ULM  =  ABS(REAL(UL(J  1 )))  +  ABS(  AIM  AG(UL(J  1 ))) 

J1=J1+N 
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IF(RN-ULM)  1,2,2 
1 RN=ULM 
2  CONTINUE 
SCL(I)=1./RN 
5  CONTINUE 
NM1=N-1 
K2  =  0 

DO  17  K=1,NM1 
BIG  =  0. 

DO  11  I  =  K,N 
IP=IPS(I) 

IPK=IP  +  K2 

SIZE  =  ( ABS(RE  AL(UL(IPK)))  +  ABS(AIM  AG(UL(IPK))))*SCL{IP) 

IF(SIZE-BIG)  11,11,10 

10  BIG  =  SIZE 
IPV  =  I 

11  CONTINUE 
IF(IPV-K)  14,15,14 

14  J  =  IPS{K) 

IPS(K)  =  1PS(1PV) 

IPS(IPV)  =  J 

15  KPP  =  IPS(K)  +  K2 
PIVOT  =  UL(KPP) 

KP1=K+1 

DO  16  I  =  KP1,N 
KP=KPP 
IP  =  IPS(1)  +  K2 
EM  =  -UL(IP)/PIVOT 
18  UL(IP)  =  -EM 
DO  16  J  =  KP1,N 
IP  =  IP  +  N 
KP  =  KP^N 

UL(1P)  =  UL(IP)  +  EM*UL(KP) 

16  CONTINUE 
K2  =  K2  +  N 

17  CONTINUE 
RETURN 
END 

SUBROUTINE  ZCURVED(NP,RO,PH,DEL,PA,NG,XG,AG,A,Z,ZMAX) 

C 

C  IMPEDANCE  ELEMENTS  FOR  CURVED  BASIS  FUNCTIONS. 

C  SPECIFICALLY  DERIVED  FOR  A  CIRCULAR  LOOP  -  NOT  A  GENERAL  CASE. 
C 

COMPLEX  CEXP,Z(50000),CON,CMPLX,SUMA,SUMB 
COMPLEX  U0,SUM  1  ,SUM2,SUM3,SUM4,EXP,ZT(500l 
DIMENSION  DEL(500),PH(500),XG(128),AG(128),PA(500) 

C  OPEN(2,FILE  =  'ZCURV  DAT  ) 

ETA  =  377. 

ZM.AX  =  0. 

PI  =  3. 14159 
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BK=2.*PI 

BK2=BK**2 

U0=(0.,0.) 

CON  =  (0. ,  1  .)*BK*ETA/(4.  *PI)'*'R0**2 
NT  =  NP-2 
C 

C  COMPUTE  Z(1,LQ)  =  ZT(LQ) 

C 

KQ=1 

Pl=DEL(KQ)/2. 

P2=PA(KQ) 

P3=DEL(KQ+l)/2. 

P4=PA(KQ+1) 

C 

C  DO  THE  L  LOOP 
C 

DO  600  LQ=1,NT 
PPl=DEL(LQ)/2. 

PP2  =  PA(LQ) 

PP3  =  DEL(LQ+l)/2. 

PP4  =  PA(LQ+1) 

C 

C  DO  THE  PHI  INTEGRATION 
C  FIRST  PART  FROM  PHI(K)  TO  PHI(K+ 1) 

C  PHI  PRIMED  INTEGRATION  FOR  THE  POSITIVE  SLOPE  OF  LQ 
C 

SUMA=U0 

PHA-PA(KQ) 

DO  100  1=1. NG 
PHI  =  P1*XG(1)  +  P2 
TK  =  (PHI-PH(KQ))/DEL(KQ) 

TKP=L/DEL(KQ)/R0 
SUM]=U0 
DO  90  J  =  1,NG 
PHIP=PP1*XG(J)  +  PP2 
TL=(PHIP-PH(LQ))/DEL(LQ) 

TLP=1./DEL(LQ)/R0 
CC  =  COS(PHI-PHIP) 

C 

C  COMPUTE  THE  MAGNITUDE  OF  R.  NOTE  THAT  R  IS  COMPUTED  FROM  THE 
C  WIRE  AXIS  TO  THE  SURFACE  OF  THE  WIRE 
C 

RR  =  R0*SQRT(4.  *(SIN((PHI-PHlP)/2.))**2  +  (A/R0)**2) 

EXP = CEXP(CM  PLX{0.  .-BK*RR))/RR 

SUM  1  =  SUM  1  +  AG(J)'*'EXP*(TK*TL*CC-TKP’»TLP/BK2) 

90  CONTINUE 

SUM1=SUM1*PP1 

C 

C  PHI  PRIMED  INTEGRATION  FOR  THE  NEGATIVE  SLOPE  OF  LQ 
C 
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SUM2=U0 
DO  80  J  =  1,NG 
PHIP=PP3*XG(J)  +  PP4 
TL  =  1 .  -(PHIP-PH(LQ  + 1  ))/DEL(LQ  + 1 ) 

TLP=-1./DEULQ+1)/R0 
CC  =  COS(PHl-PHIP) 

C 

C  COMPUTE  THE  MAGNITUDE  OF  R.  NOTE  THAT  R  IS  COMPUTED  FROM  THE 
C  WIRE  AXIS  TO  THE  SURFACE  OF  THE  WIRE 
C 

RR  =  R0*SQRT(4.  *(SIN((PHI-PHIP)/2  .))**2  +  ( A/R0)**2) 

EXP = CEXP(CMPLX(0.  .-BK*RJR))/RR 

SUM2  =  SUM2  +  AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 

80  CONTINUE 

SUM2  =  SUM2*PP3 

SUM  A  =  SUM  A  +  (SUM  1  +  SUM2)*AG(I) 

100  CONTINUE 

SUMA  =  SUMA*P1 
C 

C  SECOND  PART  FROM  PHI(K+n  TO  PHI(K  +  2) 

C  PHI  PRIMED  INTEGRATION  FOR  THE  POSITIVE  SLOPE  OF  LQ 
C 

SUMB=U0 
PHA  =  PA(K0+1) 

DO  101  1=1, NG 
PHI  =  P3*XG(1)  +  P4 
TK=  l.-(PHl-PH(KQ-  1))'DEL(KQ*  1) 

TKP  =  -l./DEL(KQ-rl)/R0 
SUM3  =  U0 
DO  91  J=1.NG 
PH1P  =  PP1*XG(J)  +  PP2 
TL  =  (PHIP-PH(LQ))  DEL(LQ) 

TLP=1./DEL(LQ)/R0 

CC=COS(PHI-PHIP) 

C 

C  COMPUTE  THE  MAGNITUDE  OF  R.  NOTE  THAT  R  IS  COMPUTED  FROM  THE 
C  WIRE  AXIS  TO  THE  SURFACE  OF  THE  WIRE 
C 

RR  =  R0*SQRT(4.  *(SlN((PHI-PHlP)/2.  ))**2  +  ( A/R0)**2) 

EXP  =  CEXP(CMPLX(0.  ,-BK*RR))/RR 

SUM3  =  SUM3  +  AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 

91  CONTINUE 

SUM3  =  SUM3*PP1 
C 

C  PHI  PRIMED  INTEGRATION  FOR  THE  NEGATIVE  SLOPE  OF  LQ 
C 

SUM4  =  U0 
DO  81  J=1.NG 
PHIP  =  PP3*XG(J)  +  PP4 
TL=  1.-(PH1P-PH(LQ+  1))/DEL(LQ  +  1) 
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493 

TLP=-1./DEL(LQ+1)/R0 

494 

CC=COS(PHI-PHIP) 

495 

C 

496 

C  COMPUTE  THE  MAGNITUDE  OF  R.  NOTE  THAT  R  IS  COMPUTED  FROM  THE 

497 

C  WIRE  AXIS  TO  THE  SURFACE  OF  THE  WIRE 

498 

C 

499 

RR=R0*SQRT(4.*(SIN((PHI-PHIP)/2.))*‘^+(A/R0)**2) 

500 

EXP = CEXP(CMPLX(0.  .-BK*RR))/RR 

501 

SUM4 = SUM4  +  AG(J)*EXP*(TK-TL*CC-TKP-TLP/BK2) 

502 

81 

CONTINUE 

503 

SUM4=SUM4*PP3 

504 

SUMB = SUMB  +  (SUM3  +  SUM4)*AG(I) 

505 

101  CONTINUE 

506 

SUMB=SUMB*P3 

507 

ZT(LQ) = CON  *(SUM  A  +  SUMB) 

508 

ZM  AX  =  AMAXKZM  AX.CABS(ZT(LQ))) 

509 

600  CONTINUE 

510 

ZT(NT)  =  ZT(2) 

511 

C 

512 

C  FILL  THE  ENTIRE  Z  MATRIX  USING  SYMMETRY  PROPERTIES 

513 

C  [Z]  IS  A  SYMMETRICAL  TOEPLITZ  MATRIX 

514 

C  ROW  INDEX,  I;  COL  INDEX,  J 

515 

C 

516 

DO  10  1=1,NT 

517 

DO  10  J  =  1,NT 

5l8 

K  =  (I-1)*NT+J 

519 

Z(K)=U0 

520 

U  =  1ABS(I-J) 

521 

IF(IJ.GT.NT)  GO  TO  10 

522 

IJ1=IJ+1 

523 

Z(K)  =  ZT(1J1) 

524 

10 

CONTINUE 

525 

RETURN 

526 

END 

527 

SUBROUTINE  CURVEW(NP,RO,PH,DEL,PA.NG.XG.AG.THR.PHR,R) 

528 

c 

529 

c 

PLANE  WAVE  EXCITATION  VECTOR  ELEMENTS  FOR  A  LOOP  USING  CURVED 

530 

c 

BASIS  FUNCTIONS.  INCIDENCE  DIRECTION  IS  (THR.PHR).  THE  WIRE  LIES 

531 

c 

IN  THE  X-Y  PLANE  (Z  =  0), 

532 

c 

533 

COMPLEX  U0,R(500),SUM  1  ,SUM2,SUM3,SUM4,CEXP.FF 

534 

COMPLEX  GG.CMPLX.SUMT.SUMP 

535 

DIMENSION  PH(500),DEL(500).PA(500),XG(128),AG(128) 

536 

MT=NP-2 

537 

U0  =  (0.,0.) 

538 

PI  =  3. 14159 

539 

BK  =  2.*P1 

540 

ST=SIN(THR) 

541 

CT  =  COS(THR) 

542 

DO  50  IP=1.MT 
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SUM1=U0 
SUM2=U0 
SUM3=U0 
SUM4  =  U0 
Pl  =  DEL(IP)/2. 

P2  =  PA(IP) 

P3=DEL(IP  +  l)/2. 

P4  =  PA(IP+1) 

DO  20  I=1,NG 
PHI=P1*'XG(I)  +  P2 
CC=COS(PHR-PHI) 

SS  =  SIN(PHR-PHI)*CT 

FF  =  AG(I)*(PHI-PH(IP))/DEL(IP)*CEXP(CMPLX(0.  ,BK*R0*ST*CC)) 

SUM1  =  SUM1  +  CC*FF 

SUM2=SUM2  +  SS*FF 

PHI  =  P3*XG(I)  +  P4 

CC  =  COS(PHR-PHI) 

SS  =  SIN(PHR-PHI)*CT 

GG  =  AG(I)*(  1 .  -(PHI-PH(1P  + 1  ))/DEL(IP  + 1  ))*CEXP(CMPLX(0. , 

*  BK*R0*ST*CC)) 

SUM3  =  SUM3-rCC*GG 
SUM4  =  SUM4  +  SS*GG 
20  CONTINUE 

SUMP  =  SUM  1  1  +SUM3*P3 

SUMT=SUM2'*P1  +SUM4*P3 
C 

C  R-WIRE-THETA  IN  R(IP)  AND  R-WIRE-PHI  IN  R(IP  +  MT) 

C 

R(IP)  =  SUMT*R0 
R(1P  +  MT)  =  SUMP*R0 
50  CONTINUE 
RETURN 
END 
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C  MAIN  PROGRAM  *LOOP.FOR* 

C  PLANE  WAVE  SCATTERING  FROM  A  CIRCULAR  LOOP  IN  THE  Z  PLANE. 

C  •***  RADAR  CROSS  SECTION  CALCULATION  ***** 

C 

COMPLEX  Z(15000).B(500).C(500).BP(500).CP(500),U 

COMPLEX  ET,EP.UC,RW(500),U0 

DIMENSION  ECP(400),IPS(500),ANG(400).EDP(400) 

DIMENSION  ZH(200),XT(128),AT(128).XH(200),YH(200) 

DIMENSION  ECV(400),EXV(400),PHC(400),PHX(400) 

DATA  PI/3.14159265/ 

DATA  IPRINT/1/ 

C 

C  READ  INPUT  AND  PROGRAM  CONTROL  PARAMETERS 
C 

OPEN(l  ,FILE=  ’PARAMLST.DAT’) 

READ(l,*)ANGLE,DT,START,STOP,AW,RB,SEG.HARR,GCONST.XLOC.XIPOL 
LOC  =  INT(XLOC) 

IPOL  =  INT(XlPOL) 

CLOSE(l) 

C 

C  READ  GAUSSIAN  CONSTANTS 
C 

OPEN(2.FILE= ’OUTGLEG') 

READ(2.*)  NT 
DO  2  1  =  1, NT 
RE.\D(2,*)  XT(I),AT(I) 

2  CONTINUE 
RAD  =  PI/180. 

ECX  =  0. 

BK  =  2.*PI 
ETA  =  377. 

U  =  (0.,1.) 

U0  =  (0..0.) 

UC  =  -U*ETA*BK'4./PI 
CONST=  16.*PI**3 
NT2  =  NT/2 
C 

C  OUTLOOP  IS  THE  FILE  THAT  THE  SCATTERED  FIELD  D.ATA  IS  WRITTEN  TO 
C  ISTOUT  IS  THE  FILE  TH.^T  THE  LOOP  CURRENT  DATA  IS  WRITTEN  TO 
C 

OPEN(8,FILE  = -OUTLOOP. DAT') 

OPEN(7,FILE=  ’ISTOUT.DAT') 

DPH1  =  SEG/RB 

NP  =  INT(2.*P1  DPHI)+  1 

DPHI  =  360./FLOAT(NP) 

W'RITE(6.*)  ’DPHI.NP=’,DPHI,NP 
C 

C  GENERATE  THE  LOOP  POINTS.  MULTIPLY  ALL  QUANTITIES  BY  BK  (  =  2*P1) 

C 

C 
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C  CHOOSE  THE  NUMBER  OF  POINTS  BASED  ON  THE  VALUE  OF  SEG 
C 

AK  =  AW*BK 
DO  10  I=1,NP+1 
PP=FLOAT(l-]  )*DPHI*RAD 
XH{I)  =  RB*COS(PP)*BK 
YH(I)  =  RB*S1N(PP)*BK 
ZH(I)=0. 

10  CONTINUE 
NP=NP+2 
C 

C  OVERLAP  THE  ENDS  SO  THAT  THE  CURRENT  IS  CONTINUOUS 
C 

XH(NP)=XH(2) 

YH(NP)=YH{2) 

ZH(NP)  =  ZH(2) 

DO  52  I  =  1.NP 
YHB  =  YH(1)/BK 
XHB  =  XH(I)'BK 
DEL  =  0. 

IPd.NE.DTHEN 
DXX=XHB-XH(I-1)/BK 
DYY  =  YHB-YH(I-1)/BK 
DEL  =  SQRT(DXX**2  -r  DYV*2) 

ENDIF 

52  CONTINUE 
C 

C  DEFINE  DIMENSIONS  OF  THE  IMPEDANCE  MATRIX  BLOCKS 
C 

MT  =  NP-2 

WR1TE(6.*)  'MT^dMT 
C 

C  COMPUTE  IMPEDANCE  MATRIX  ELEMENTS 
C 

CALL  ZMA'TVSVVd.LNP.RB.XH.YH.ZH.NT.XT.AT.AK.Z) 

WRITE(6,*)  'Z  COMPUTED' 

1F(IPRINT.EQ.0)THEN 
DO  1010  1=1. MT 
CZ  =  CABSIZ(I)) 

AZ  =  ATAN2(A1MAG(Z(1)).REAL(Z(1))+  l.E-8)/RAD 
WR1TE(6,*)  '1,Z=  '.l.ZdKCZ.CZ/ZMAX.AZ 
1010  CONTINUE 
ENDIF 
C 

C  PERFORM  LU  DECOMPOSITION 
C 

CALL  DECOMP(MT.lPS,Z) 

WRITE(6,*)  'Z  DECOMPOSED' 

C 

C  BEGIN  FIELD  CALCULATIONS  PHI  FOR  PATTERN  CUT  (DEGREES)  =  PHID 
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C 

PHID=0. 

PHR0=PH1D*RAD 
IT = INT((STOP-START)/DT)  + 1 
WRITE(7,*)  rr.MT.o.o 
DO  500  1  =  1, IT 

THETA = FLO  AT(I  - 1  )*DT + START 

THR=THETA*RAD 

PHR=PHR0 

IFCTHETA.LT.  180.)  GO  TO  99 
THR  =  (360.  -THETA)*RAD 
PHR=PHR0+PI 
99  CONTINUE 
ET=U0 
EP=U0 
C 

C  COMPUTE  THE  EXCITATION  VECTOR 
C 

CALL  PLANEW(NP.XH.YH,ZH.THR.PHR,RW) 

IF  (LOC  .EQ.  0)  THEN 
C 

C  CIRCULAR  POLARIZATION  IF  LOC=l  ELSE  LINEAR 
C 

IF(IPOL.EQ.l)THEN 

C 

C  THETA  POLARIZED  INCIDENT  WAVE  (IPOL=  1) 

C 

DO  101  L=l,MT 

101  B(L)  =  RW(L) 

ELSE 

C 

C  PHI  POLARIZED  INCIDENT  WAVE  (IPOL  =  2) 

C 

DO  102  L=1,MT 

102  B(L)  =  RW(L  +  MT) 

ENDIF 

C 

C  PERFORM  GAUSSIAN  ELIMINATION  TO  DETERMINE  (CJ 
C 

C.4LL  SOLVE(MT,lPS.Z,B,C) 

DO  210  L=1.MT 

WR1TE(7,*)  L.L*DPHLCABS(C(L)).ATAN2(REAL(C(L)). 
*  AIMAG(C(L)))/RAD 

ET = ET  +  (RW(L)/BK)*C(L) 

210  EP=EP  +  (RW(L  +  MT)/BK)*C(L) 

ELSE 

C 

C  THETA  PLOARIZED  INCIDENT  WAVE 
C 

DO  221  L=1.MT 
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221  B(L)  =  RW(L) 

C 

C  PHI  POLARIZED  INCIDENT  WAVE 
C  PHASE  SHIFT  FOR  CP  IS  PI/2. 

C 

DO  222  L=1,MT 

222  BP(L)=RW(L  +  MT)*CEXP{CMPLX(0..PI/2.)) 

CALL  SOLVE(MT,IPS.Z.B.C) 

CALL  SOLVE(MT,IPS,Z.BP.CP) 

DO  220  L=I,MT 

WRITE(7.*)  L.L*DPHI,CABS(C(L)  +  CP(L)),ATAN2(REAL(C(L)  +  CP(L)). 
*  AIMAG(C(L)  +  CP(L)))/RAD 

ET = ET + RW(L)*C(L)  +  R  W(L)*CP(L) 

220  EP  =  EP  +  RW(L  +  MT)*C{L)  +  RW^L  +  MT)*CP(L) 

ENDIF 

ET=UC*ET 

EP=UC*EP 

C 

C  E-THETA  IS  CO-POL;  E-PHI  IS  CROSS-POL 
C 

ANGfl)  =  THETA 
ECV(I)  =  CABS(ET) 

EXV(I)  =  CABS(EPl 
ECR  =  REAL(ET) 

ECI  =  AIMAG(ET) 

EXR  =  REAL(EP) 

EXI  =  AIMAG(EP) 

PHC(l)  =  ATAN2(ECI.ECR  +  l.E-20>/RAD 
PHX(I)  =  ATAN2(EXLEXR-H.E-20)'RAD 
ECX  =  AMAX1(ECX.ECV(I).EXV(I)) 

500  CONTINUE 
DO  600  1=1,IT 

ECV(I)  =  AMAX1(ECV(I).1.E-10) 

EXV(I)  =  AMAX1(EXV(I),1.E-10) 

ECP(l)  =  (ECV(I)/ECX)**2 
EDP(I)  =  »EXV(I)/ECX)**2 
ECP(1 1  =  AM  AX  1  (ECP(I ),  .0000 1 ) 

EDP(I  )  =  AM  AXKEDPd).. 00001) 

ECP(1)=  10.*ALOG10(ECP(1)) 

EDP(I)=  10.*ALOG10(EDP(I)) 

600  CONTINUE 

SIGM  A  =  (ECX**2)*C  ABS(UC).'(2.*BK) 

S1GDB=10.*ALOG10(SIGMA) 

WR1TE(6.*)  'SIGMA,  IN  DB  = '.MGMA.SIGDB 
DO  9000  L=1.IT 
WRITE(8,*)  ANG(L).ECV(L) 

9000  CONTINUE 
900  CONTINUE 
STOP 
END 


SUBROUTINE  SOLVE(N,IPS.UL.B,X) 
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C 

C  SUBROUTINE  TO  SOLVE  SYSTEM  OF  EQUATIONS  WITH  COMPLEX 
COEFFICIENTS. 

C  CALL  ’DECOMP’  FIRST.  (FROM  MAUTZ  AND  HARRINGTON) 

C 

COMPLEX  UL(50000),B(500),X(500).SUM 
DIMENSION  IPS(500) 

NP1=N  +  1 
IP=IPS(1) 

X{1)=B(IP) 

DO  2  1=2, N 
IP=IPS(I) 

IPB  =  IP 
IM  1  =  1-1 
SUM  =  (0.,0.) 

DO  1  J=1,IM1 

SUM  =  SUM-hUL(IP)*X(J) 

1  1P=IP  +  N 

2  X(I)  =  B(IPB)-SUM 
K2  =  N*(N-1) 

IP=IPS(N)  +  K2 
X(N)  =  X(N)/UL(IP) 

DO  4  IBACK  =  2,N 
1  =  NP1-1BACK 

K2  =  K2-N 
IPl  =  IPS(I)-hK2 
lPl  =  l+l 
SUM  =  (0..0.) 

IP  =  IPI 

DO  3  J  =  IP1.N 
IP  =  IP-i-N 

3  SUM  =  SUM  +  UL(IP)*X(J) 

4  X(1)  =  {X{I)-SUM)/UL(1P1) 

RETURN 

END 

SUBROUTINE  DECOMP(N,IPS.UL) 

C 

C  SUBROUTINE  TO  DECOMPOSE  SYSTEM  OF  EQUATIONS. 

C  FROM  MAUTZ  AND  HARRINGTON. 

C 

COMPLEX  UL(50000),P1VOT.EM 
DIMENSION  SCL(500),IPS(500) 

DO  5  I=1,N 
IPS(I)  =  I 
RN=0. 

J1  =  I 

DO  2  J  =  1.N 

ULM  =  ABS(RE  AL(UL(J  I )))  +  ABS(  AIM  AG(UL(J  1 ))) 

Jl=Jl-(-N 
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IF(RN-ULM)  1,2,2 
1 RN=ULM 
2  CONTINUE 
SCL(1)=1./RN 
5  CONTINUE 
NM1=N-1 
K2=0 

DO  17  K=1,NM1 
BIG=0. 

DO  11  I=K,N 
IP=IPS(I) 

IPK=IP  +  K2 

SIZE  =  (ABS(REAL(UL(IPK)))  +  ABS(  AIM  AG(UL(1PK))))*SCL(1P) 

IF(SIZE-BIG)  11,11,10 

10  BIG  =  SIZE 
IPV  =  I 

11  CONTINUE 
IF(IPV-K)  14,15,14 

14  J  =  1PS(K) 

IPS(K)  =  IPS(IPV) 

IPS(IPV)  =  J 

15  KPP  =  IPS(K)  +  K2 
PIVOT  =  UL(KPP) 

KP1=K+1 

DO  16  1  =  KP1,N 
KP  =  KPP 
1P  =  1PS(I)  +  K2 
EM  =  -UL('P)/PIVOT 
18  UU1P)  =  -EM 
DO  16  J  =  KP1.N 
IP  =  1P  +  N 
KP  =  KP  +  N 

UL(IP)  =  UL(IP)  +  EM*UL(KP) 

16  CONTINUE 
K2  =  K2  +  N 

17  CONTINUE 
RETURN 
END 

SUBROUTINEMATW\\(NWlRES,NWl,NW2,RB.XH.YH,ZH,NT,XT,AT,AK,ZZl 

C 

C  IMPEDANCE  ELEMENTS  FOR  LINEAR  BASIS  FUNCTIONS. 

C 

COMPLEX  CEXP,Z(200).ZZ(15000),CON,CMPLX,EXP 

COMPLEX  U0,SUM.SUM1,SUM2.SUM3.SUM4 

DIMENSION  ZH(200),XT(128).AT(128),XH(200),YH(200).UU(200) 

DIMENSION  D1(200),S1(200),C1(200),ZS1(200) 

DIMENSION  XS1(200).YS1(200) 

DIMENSION  CU(200).SU(200) 

INTEGER  NT.NWIRES.NW1(4).NW'2(4),NS(4) 

PI  =  3. 1415926 


59 


60 


351 

352 

353 

354 

355 

356 

357 

358 

359 

360 

361 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 


X1  =  XS1(1)+T*CU(I) 

YI  =  YS1(I)+T*SU(I) 

ZI  =  ZS1(I) 

DO  100  L=1,NT 
TP  =  T2*XT(L) 

TJ  =  .5+TP/D1(J) 

XJ=XS1(J)+TP*CU(J) 

YJ  =  YS1(J)  +  TP*SU(J) 

ZJ  =  ZS1(J) 

RP = SQRT((XI-XJ)**2  +  (YI-YJ)*-»2  +  A**2) 

EXP = CEXP(CMPLX(0.  .-RP))/RP 

SUM =SUM  +  AT(L)*AT(K)*EXP*(TI*TJ*CC-TIP*TJP) 

100  CONTINUE 

SUMl  =SUM-Tl*CON-^ 

C 

C  DOING  12 
C 

J=JQ+  1 

CC  =  COS(UU(I)-UU(J)) 

TIP=1./D1(D 
TJP  =  -1./D1(J) 

Tl=DI(I)/2. 

T2  =  Dl(J)/2. 

SUM  =  U0 
DO  101  K=1,NT 
T=T1*XT(K) 

TI  =  .5  +  T/D1(1) 

XI  =  XSI(I)  +  T*CL'(f) 

YI  =  YSl(l)-rT*SU(l) 

Z1  =  ZS1(I) 

DO  101  L=  l.NT 
TP  =  T2*XT(Li 
TJ  =  .5-TP/D1(J) 

XJ  =  XS1(J)  +  TP*CU(J) 

YJ  =  YS1(J)  +  TP*SU(J) 

ZJ  =  ZS1(J) 

RP  =  SQRT((X1-XJ)**2  +  (Y1-YJ)**2  +  A**2) 

EXP  =  CEXP(CMPLX(0..-RPl)/RP 
SUM=SUM+AT(L)*AT(K)*EXP*(T1'*‘TJ*CC-T1P*TJP) 

101  CONTINUE 

SUM2  =  SUM  *T  1  *CON  *72 
C 

C  DOING  13 
C 

I  =  1P+1 
J  =  JQ 

CC  =  COS(UU(I)-UU(J)) 

TIP  =  -1./D1(1) 

TJP=1./D1(J) 

Tl=Dl(l)/2, 
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Z(JQ)  =  (SUM  1  +  SUM2  +  SUM3  +  SUM4) 

600  CONTINUE 
Z(NTR1A)  =  Z(2) 

C 

C  FILL  THE  ENTIRE  Z  MATRIX  USING  SYMMETRY  PROPERTIES 
C  ROW  INDEX,  i;  COL  INDEX,  J 
C 

DO  12  I=1,NTRIA 
DO  12  J=1.NTRIA 
K  =  (1-1)*NTRIA  +  J 
ZZ(K)  =  U0 
U  =  1ABS(1-J) 

IF(U.GT.NTRIA)  GO  TO  12 

U1=U+1 

ZZ(K)  =  Z(1J1) 

12  CONTINUE 
CLOSE(2) 

RETURN 

END 

SUBROUTINE  PLANEW(NP,XH.YH,ZH,THR,PHR,R) 

C 

C  PLANE  WAVE  EXCITATION  VECTOR  ELEMENTS  FOR  WIRE  AND 
C  INCIDENCE  DIRECTION  IS  (THR.PHR). 

C  WIRE  LIES  IN  THE  X-V  PLANE  (Z  =  0) 

C 

COMPLEX  U0.C.R(2000).CEXP,EXP,F11,F12.SI,DI,CMPLX 
DIMENSION  ZH(500).XH(500).YH(500) 

MP2  =  NP-! 

MT2  =  NP-2 
U0  =  (0..0. ) 

CC  =  COS(THR) 

SS=SIN{THR) 

CP  =  COS(PHR) 

SP  =  SIN(PHR) 

UP  =  SS*CP 
VP=SS*SP 
DO  12  IP=1.MP2 
I1  =  IP 
1  =  IU  1 

ZS  =  .5*(ZH(I)  +  ZH(II)) 

XS=.5*(XH(I)  +  XH(II)) 

YS  =  .5*(YH(I)^  YH(II)) 

DX  =  XH(1)-XH(II) 

DY  =  VH(I)-YH(II) 

D1=SQRT(DX**2  +  DY**2) 

SU  =  DY/D1 
CU  =  DX/D1 

C  FOR  WRES  IN  THE  XY  P  ,ANE  SL  (V)=  I  AND  COS(V)  =  0 
SV=1.0 
CV  =  0.0 
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C 

C  WIRE  SEGMENT  CALCULATIONS 
C 

A  =  UP*CU  +  VP*SU 
B=UP*XS  +  VP*YS 
C  =  CMPLX(0..A) 
EXP=CEXP(CMPLX(O..B)) 

AA  =  CC*(CU*CP  +  SU*SP) 

BB  =  SU*CP-SP*CU 
PSl=Dl*A/2. 

1F(PSI.NE.0.)  GO  TO  60 
SINC=1. 

GO  TO  61 

60  SINC  =  S1N(PSI)/PSI 

61  COSP=COS(PSI) 

FIl  =  SINC*Dl*EXP/2. 

F12  =  (0.,0.) 

IF(ABS(A).LT.  1  .E-4)  GO  TO  62 
CSP  =  COSP-SINC 
IF(ABS(CSP).LT.  1  .E-4)  GO  TO  62 
FI2  =  EXP/C*CSP 

62  CONTINUE 
SI  =  FIl-t-FI2 
D1  =  F11-FI2 

C 

C  R-WIRE-THETA 
C 

IF(IP.EQ.MP2)GO  TO  10 
RaP)  =  AA-Sl 
R(1P  +  MT2)  =  BB*SI 
10  CONTINUE 
C 

C  R-WIRE-PHI 
C 

14  IFdP.EQ.l)  GO  TO  12 
R(1P-1)  =  R(IP-1)-^AA*D1 
R(IP- 1  -t-  MT2)  =  R(IP-1  -t-  MT2)  +  BB*DI 
12  CONTINUE 
210  RETURN 
END 
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C  MAIN  PROGRAM  *HARLOOP.F* 

C  PLANE  WAVE  SCATTERING  FROM  A  CIRCULAR  LOOP  IN  THE  Z  PLANE. 

C  ****  RADAR  CROSS  SECTION  CALCULATION  ***** 

C  USING  HARRINGTON'S  FORMULATION  FROM  THE  BOOK  'FIELD  COMP.  BY 
C  MM' (P.83  TO  95) 

C 

COMPLEX  Z(5000).E(250),C(250).EX.EC.ET,EP,RW(1000),UC,EPHI(500) 
COMPLEX  CPHI(500) 

DIMENSION  ECP(500),IPS(250).ANG(500),EDP(500),XT(300).AT(300) 
DIMENSION  ECV(500),EXV(500),PHC(500),PHX(500) 

DATA  PI/3.14159265/ 

DATA  IPRINT/1/,ITEST/1/ 

RAD  =  PI/180. 

ECX  =  0. 

BK=2.*PI 
ETA  =  377. 

UC  =  (0,- 1  )*ETA*BK/(4.  *PI) 

PHIRD  =  0. 

OPEN(l  ,FILE  =  'PARAMLST.DAT') 

READ(I,*J  ANGLE, DT, START, STOP. A.B.SEG.AHARR.GCONST.XLOC.XIPOL 
1P0L  =  1NT(XIP0L) 

L0C  =  INT(XL0C) 

CLOSE(l) 

NM  =  1NT(AHARR) 

CON  =  {377.*BK)**2/2./BK 
0PEN(2,FILE  =  '0UTGLEG') 

RE.AD(2,*)  NT 
DO  I  1=1. NT 
READ(2,*)  XT(I).AT(I) 

1  CONTINUE 

0PEN(8.FILE  =  '0UTH,ARR.D,AT') 

OPt.Nf7.F[LE= 'IHARROUT.D.AT') 

NR0W  =  2*NM-^  1 
WRITE(6,1300)  B.A.NROW.NT 

1300  FORM.AT(//,5X.'PL.XNE  WAVE  SCATTERING  BY  A  CIRCULAR  LOOP',/ 

*  .5X,  'USING  METHOD  IN  HARRINGTONS  MM  TEXT  BOOK',/ 

*  .5X,'LOOP  RADIUS  (WAVL)= '.F8.4,/,5X,'WIRE  R.ADIUS  (WAVL)=', 

*  F8. 4, /,5X. 'NUMBER  OF  AZIMUTHAL  MODES  (INCLUDING  ZERO)=',I3. 

*  /,5X, 'NUMBER  OF  INTEGRATION  POINTS  IN  PHI  = ',14) 

C 

C  COMPUTE  IMPEDANCE  MATRIX  ELEMENTS 
C 

W'RITE(6,*)  'CALLING  ZMAT' 

CALL  ZMATW\\'(NM.A.B.NT,XT.AT.Z) 

WR1TE(6,*)  'WIRE  IMPEDANCE  COMPUTED' 

CALL  DECOMP(NROW.IPS.Z) 

WR1TE(6,*)  'Z  DECOMPOSED' 

C 

C  BEGIN  FIELD  CALCULATIONS 
C 
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PHR0=PHIRD*RAD 
rr = INT((STOP-START)/DT)  +  1 
WRITE{7.*)  IT,NROW 
DO  500  1  =  1, IT 

THETA  =  FLOAT(I- 1  )*DT  +  START 
WRITE(6,*)  THETA  = ’.THETA 
THR  =  THETA*RAD 
PHR=PHR0 

IF(THETA.LT.180.)  GO  TO  99 
THR  =  (360.-THETA)*RAD 
PHR  =  PHR0+P1 
99  CONTINUE 
ET  =  (0..0.) 

EP  =  (0.,0.) 

CALL  PLANEW(NM,B,THR.PHR,RW) 

C  TRANSMIT  VECTOR  ELEMENTS  ARE  TRANSPOSED  FORMS  OF  RECEIVE 
VECTOR 

C  ALSO  THE  THETA  COMPONENT  GETS  A  NEGATIVE  SIGN 
lF(LOC  .EQ.  0)  THEN 
IFdPOL.EODTHEN 
C 

C  THETA  POLARIZED  INCIDENT  WAVE  (IPOL=  1) 

C 

DO  101  L=l,NROW 

101  E(NR0W-L+1)  =  RW(L) 

ELSE 

C 

C  PHI  POLARIZED  INCIDENT  WAVE  (1P0L=2; 

C 

DO  102  L=  LNROW 

102  E(NROW-L+ l)  =  RW(L  +  NROW) 

ENDIF 

WRITE(6.*)  ’CALLING  SOLVE’ 

CALL  SOLVEINROW.IPS.Z.E.C) 

WRITE(6,*)  ’RETURNED  FROM  SOLVE’ 

DO  210  L=  LNROW 
WR1TE(7.*)  C(L) 

ET  =  ET  +  RW(L)*C(L) 

210  EP  =  EP  +  RW(L  +  NROW)*C(L) 

ELSE 

C 

C  E-THETA  IS  CO-POL;  E-PHl  IS  CROSS-POL 
C 

DO  221  L=  LNROW 

221  E(NROW-L-t-l)=RW(L) 

DO  222  L=  LNROW 

222  EPHI(NROW-L-t-l)  =  RW(L-t-NROW)*CEXPlCMPLX(0.,Pl/2.)) 

CALL  SOLVE(NROW.IPS,Z.E.C) 

CALL  SOLVE(NROW,IPS,Z.EPHl.CPHI) 

DO  220  L=  LNROW 


66 


101 

102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 


WRITE(7,*)  C(L)  +  CPHI(L) 

ET = ET + R  W(L)-^C(L)  +  R  W(L)*CPHI(L) 

220  EP  =  EP  +  RW{L  +  NROW)*C(L)  +  R  W(L  +  NROW)  *CPHI(L) 

ENDIF 
EC=UC*ET 
EX=UC*EP 
ANG(I)=THETA 
ECV(I)  =  CABS(EC) 

EXV(1)  =  CABS(EX) 

ECR=REAL(EC) 

ECI  =  AIMAG(EC) 

EXR=REAL(EX) 

EXI  =  AIMAG(EX) 

PHC(l)  =  ATAN2(ECI,ECR  + 1  .E-20)/R  AD 
PHX(I)  =  ATAN2(EX1 ,  EXR  + 1 .  E-20)/RAD 
ECX  =  AMAX1(ECX,ECV(I),EXV(1)) 

500  CONTINUE 
DO  600  I=1,IT 

ECV(I)  =  AM  AX  1  (ECV(I),  1  .E- 10) 

EXV(I)  =  AMAX1(EXV(I).1.E-10) 

ECP(1)  =  {ECV(I)/ECX)**2 
EDP(I)  =  (EXV(I)/ECX)**2 
ECP(I)  =  AM  AX  1(ECP(I),  .0000 1 ) 

EDP(1)  =  AM  AXKEDPd).. 00001) 

ECP(I)=  10.*ALOG10(ECP(1)) 

EDP(1)=  10.*ALOG10(EDP(1)) 

600  CONTINUE 

SIGMA  =  (ECX**2)*4.*PI 
S1GDB=  10.*ALOG10(S1GMA) 

WRITE(6,*)  'BACKSCATTER,  IN  DB  = ',S1GMA,SIGDB 
OPEN(2.FILE=’RCS.DAT') 

WRITE(2,'*)  SIGMA, SIGDB 
CLOSE(2) 

C 

C  PRINT  FIELD  POINTS 
C 

DO  9000  L=  LIT 
WRITE(8,*j  ANG(L),ECV(L) 

5016  FORM  AT(5X.F6. 1 .3X,2(F8.4.3X.F7.  L3X.F7.2,3X)) 

9000  CONTINUE 
900  CONTINUE 
9998  STOP 
END 

SUBROUTINE  SOLVE(N,IPS,UL.B.X) 

COMPLEX  UL(5000).B(250).X(250),SUM 
DIMENSION  IPS(250) 

NP1=N+1 
IP  =  IPS{1) 

X(1)  =  B(IP) 

DO  2  1  =  2.N 
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IP  =  IPS(1) 

IPB=IP 
IM1  =  M 
SUM  =  (0.,0.) 

DO  1  J=1,IM1 

SUM  =  SUM+UL(1P)*X(J) 

1  IP=IP+N 

2  X(I)=BaPB)-SUM 
K2  =  N*(N-1) 

IP=1PS(N)  +  K2 
X(N)=X(N)/UL(IP) 

DO  4  IBACK  =  2,N 
1  =  NP1-IBACK 

K2  =  K2-N 
IPI  =  IPS(I)  +  K2 
IP1=1+1 
SUM  =  {0.,0.) 

IP  =  IPI 

DO  3  J  =  IP1.N 
IP  =  1P  +  N 

3  SUM  =  SUM+UL(IP)*X(J) 

4  X(1)  =  (X(I)-SUM)/UL(1P1) 

RETURN 

END 

SUBROUTINE  DECOMP(N.lPS,UL) 

COMPLEX  UL(5000),P1VOT,EM 
DIMENSION  SCL(250),IPS(250) 

DOS  I=I,N 
IFS(I)  =  I 
RN  =  0. 

J1=I 

DO  2  J=1,N 

ULM  =  ABS(REAL(UL(J  1 )))  +  ABS(  AIM  AG(UL(J  1 ))) 

J1=J1+N 
IF(RN-ULM)  1,2.2 
1 RN=ULM 
2  CONTINUE 
SCL(I)=1./RN 

5  CONTINUE 
NM1=N-I 
K2  =  0 

DO  17  K  =  1,NM1 
B1G  =  0. 

DO  11  1  =  K,N 
IP  =  IPS(I) 

IPK  =  IP  +  K2 

SIZE  =  (ABS(REAL(UL(IPK)))  +  ABS(AIMAG(UL(1PK))))*SCL(1P) 
IF(SIZE-BIG)  11,11,10 
10  BIG  =  SIZE 
IPV  =  I 
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1 1  CONTINUE 
IFaPV-K)  14,15,14 

14  J  =  IPS(K) 

IPS(K)=IPS(1PV) 

IPS(IPV)  =  J 

15  KPP=IPS(K)  +  K2 
PIVOT=UL(KPP) 

KP1=K+1 

DO  16  1=KP1,N 
KP=KPP 
IP=IPS(I)  +  K2 
EM  =  -UL(IP)/PIVOT 
18  UL(IP)  =  -EM 
DO  16  J  =  KP1,N 
IP  =  IP  +  N 
KP  =  KP  +  N 

UL{IP)  =  UL(IP)  +  EM  -UUKE) 

16  CONTINUE 
K2  =  K2  +  N 

17  CONTINUE 
RETURN 
END 

SUBROUTINE  ZMATW'Vk'(NM.A.B,NT,XT.AT,Z) 

C 

C  ***  MODS  FOR  LOOP  -  USING  HARRINGTON'S  TEXT  BOOK  EQUATIONS  AS  A 
C  CHECK  FOR  OTHER  METHODS.  NM  IS  THE  NUMBER  OF  AZIMUTHAL  MODES 
USED 
C 

COMPLEX  Z(5000).CON,CMPLX.FK 
DIMENSION  XT(300)..ATi300) 

PI  =  3. 1415926 
BK  =  2  *PI 

CON  =  CMPLXf0,.PI*377.*BK*B) 

NROW  =  2.*NM^  1 
DO  10  I  =  -NM.NM 
DO  10  J  =  -NM.NM 
IJ=J  +  NM-  1  ^(I-rNM)*NRO\V 
Z(1J)  =  (0.,0,) 

10  CONTINUE 

'  ONLY  DIAGONAL  ELEMENTS  ARE  NONZERO.  ALTHOUGH  SYMMETRY  EXISTS 
.sETANTEN 

•:  Z(-N.-N)  AND  Z(N.N)  IT  IS  NOT  BEING  USED. 

DO  20  I  =  -NM,NM 
J  =  I  +  NM  1  +  (I  +  NM)*NROW 
IP=I+I 
1M  =  I-1 

Z(J)  =  (.5*FK(1M.B,A.NT.XT,AT)+.5*'FK(IP,B.A,NT,XT.AT)- 
*  (I  BK'B)**2»FK(I.B,A,NT.XT.AT))*CON 
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20  CONTINUE 
RETURN 
END 

SUBROUTINE  PLANEW(N.B.THR.PHR,R) 

C 

C  PLANE  WAVE  RECEIVE  VECTOR  ELEMENTS  FOR  WIRE  USING  THE 
C  FORMULATION  FROM  HARRINGTON’S  BOOK.  N  IS  THE  NUMBER  OF 
C  AZIMUTHAL  MODES.  NOTE  THAT  B(N)  =  R(-N)  (B  IS  EXCITATION  AND 
C  R  IS  RECEIVE). 

C 

COMPLEX  R(1000),CEXP,EXP.CMPLX 
PI  =  3. 14159 
BK  =  2.*P1 
CT=COSrrHR) 

ST  =  SIN(THR) 

RR  =  BK*B*ST 
NROW  =  2*N  +  l 

C  DO  THETA  RECEIVE  COMPONENTS  FIRST 
DO  10  I  =  -N,N 
IP  =  I+1 
IM  =  M 

EXP  =  CEXP(CMPLX(0.,I*PHR)) 

Rf  I  -r  N  +  1 )  =  -PI*B*(0. ,  1 .  )**1*EXP*(BESSJ(1P.RR)  +  BESSJ(IM  .RR))*CT 
10  CONTINUE 

C  NOW  DO  PHI  RECEIVE  COMPONENTS 
DO  20  I  =  -N.N 
1P  =  I+1 
IM  =  I-I 

EXP  =  CEXP(CMPLX(O..I-*PHR)) 

R(I  +  N  +  1  +  NROW)  =  PI*B*(0.,1. )•*(!+  1)*EXP*(BESSJ(1P,RR) 

•  -BESSJdM.RR)) 

20  CONTINUE 
RETURN 
ENC 

COMPLEX  FUNCTION  FKfN.B.A.NT.XT.AT) 

COMPLEX  SUM.EXP1,EXP2.CEXP.CMPLX 
DIMENSION  XT(300).AT(300) 

PI  =  3. 14159 
BK  =  2.*PI 
CC=1./BK 
Pl=(2.*Pl-0.)/2. 

P2  =  (2.*Pl+0.)/2. 

SUM  =  (0.,0.) 

DO  10  1=1, NT 
PHI  =  P1*XT(1)  +  P2 

RR  =  SQRT(4.*SIN(PHI/2.)**2+(A/Br*2) 

EXPl  =CEXP(CMPLX(0.,-BK*B*RR)) 

EXP2  =  CEXP(CMPLX(0..-N*PHI)) 

SUM  =  SUM  +  AT(1)*EXP1  *EXP2/RR 
10  CONTINUE 
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FK  =  SUM*P1*CC 

RETURN 

END 

FUNCTION  BESSJ(NN,X) 

C  RETURNS  THE  BESSEL  FUNCTION  B  OF  ORDER  N  (>  1)  AND  REAL 
C  ARGUMENT  X. 

PARAMETER  (IACC  =  40,BIGNO=  I.EIO.BIGNI  =  1  .E-10) 

IF(NN.LT.O)  N=-NN 
IF(NN.GE.O)  N  =  NN 
KC  =  3 

IF(N.EQ.O)  KC=1 
IF(N.EQ.1)KC=2 
GO  TO  (1,2, 3), KC 

1  BESSJ=BESSJO(X) 

GO  TO  4 

2  BESSJ  =  BESSJ1(X) 

GO  TO  4 

3  BESSJ  =  0. 

IF(ABS(X).LT.l.E-5)  GO  TO  4 
TOX  =  2./X 

IF(X.GT.FLOAT(N)j  THEN 
BJM  =  BESSJ0(X) 

BJ  =  BESSJ1{X) 

DO  11  J  =  1,N-1 
BJP  =  J*TOX*BJ-BJM 
BJM  =  BJ 
BJ  =  BJP 

1 1  CONTINUE 
BESSJ  =  BJ 

ELSE 

M  =  2*((N  -^INT(SQRT(FLOAT(lACC*N))))/2) 

BESSJ  =  0. 

JSUM=0. 

SUM  =  0. 

BJP  =  0. 

BJ  =  1. 

DO  12  J  =  M.1.-1 
BJM=JTOX*BJ-BJP 
BJP  =  BJ 
BJ  =  BJM 

IF(ABS(BJ).GT.BIGNO)  THEN 
BJ  =  BJ*BIGNI 
BJP  =  BJP*BIGNI 
BESSJ  =  BESSJ*BIGNI 
SUM  =  SUM-»BIGN1 
ENDIF 

IF(JSUM.NE.O)  SUM  =  SUM  +  BJ 

JSUM=1-JSUM 

IF(J.EQ.N)  BESSJ  =  BJP 

12  CONTINUE 
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SUM=2.*SUM-BJ 

BESSJ=BESSJ/SUM 

ENDIF 

4  CONTINUE 

IF(NN.LT.O)  BESSJ  =  (-1 , )**N*BESSJ 

RETURN 

END 

FUNCTION  BESSJO(X) 

C 

C  BESSEL  FUNCTION  OF  0  ORDER,  REAL  ARGUMENT  X 
C  (SEE  ’NUMERICAL  RECIPES’.  P.172) 

C 

REAL*8  Y,P1,P2,P3,P4,P5,Q1,Q2.Q3.Q4,Q5.R1.R2,R3.R4,R5,R6. 

*  S1,S2,S3,S4.S5.S6 

DATA  P1.P2,P3,P4.P5/1.D0.-.109862827D-2..2734510407D-4, 

*  -.2073370639D-5. . 209388721  lD-6/ 

DATA  Q1.Q2.Q3,Q4.Q5/-.1562499995D-1,.1430488765D-3. 

*  -.6911147651D-5,.7621095161D-6,-.934945152D-7/ 

DATA  R1,R2.R3,R4,R5.R6/57568490574. DO, -13362590354. DO, 

*  651619640.7D0.-11214424.18D0.77392.33017D0,-184.9052456D0/ 
DATA  SI, S2,S3,S4,S5,S6/57568490411. DO,  1029532985. DO, 

*  9494680. 7 1 8D0,59272.64853D0,267. 85327 12D0. 1  .DO/ 
IF(ABS(X).LT.8.)THEN 

Y  =  X**2 

BESSJ0  =  (R 1  Y*(R2  Y*(R3  +  Y*(R4  +  Y'*(R5  +  Y*R6)))))/ 

*  (S 1  +  Y*(S2  +  Y*(S3  +  Y*(S4  +  Y*(S5  -h  Y*S6))))) 

ELSE 

AX=  ABS(X) 

Z=8./AX 

Y  =  Z**2 

XX  =  AX-. 785398164 

BESSJ0  =  SQRT(. 636619772  AX)*(COS(XX)nPl+Y*(P2-^Y*(P3  + 

*  Y*(P4-t-Y*P5))))-Z*SlN(XX)nQl  +  Y'‘(Q2 -r- Y*(Q3 

*  Y»(Q4  +  Y*Q5))))) 

ENDIF 

RETURN 

END 

FUNCTION  BESSJ  1(X) 

C 

C  BESSEL  FUNCTION  B  OF  ORDER  1.  REAL  ARGUMENT  X 
C  (SEE  ’NUMERICAL  RECIPES'.  P.173) 

C 

REAL*8  Y,P1.P2,P3,P4,P5.Q1.Q2.03,Q4,Q5,R1.R2,R3,R4,R5,R6, 

*  S1,S2.S3,S4.S5,S6 

DATA  P1,P2,P3.P4,P5/1.D0..183105D-2,-.35I6396496D-4, 

*  .2457520I74D-5,-.20337019D-6/ 

DATA  Q1,Q2,Q3,Q4,Q5/.04687499995D0.-.2002690873D-3. 

*  .8449199096D-5,-.99228987D-6..  105787412D-6/ 

DATA  R1,R2,R3,R4,R5,R6/723626I4232.D0,-7895059235.D0. 

*  242396853. 1D0.-2972611. 43900,15704. 4826D0,-30.16036606D0; 
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401  DATA  S1,S2,S3,S4,S5,S6/144725228442.D0,2300535178.D0, 

402  18583304. 74D0,99447.43394D0, 376.999139700,1. DO/ 

403  IF( ABS(X) .  LT.  8 . )  THEN 

404  Y  =  X**2 

405  BESSJ1  =  X*(R1  +Y*(R2  +  Y*(R3  +  Y*(R4  + Y*(R5+Y*R6)))))/ 

406  *  (S1+Y*(S2  +  Y*(S3  +  Y*(S4  +  Y*(S5  +  Y*S6))))) 

407  ELSE 

408  AX  =  ABS(X) 

409  Z  =  8./AX 

410  Y  =  Z**2 

411  XX  =  AX-2. 356 194491 

412  BESSJl  =SQRT(.636619772/AX)*(COS(XX)*(Pl  -l■Y*(P2-l■  Y*(P3-i- 

413  ♦  Y*(P4-t-Y*P5))))-Z*SIN(XX)*(QH-Y*(Q2-i-Y*(Q3-i- 

414  *  Y*(Q4-f-Y*Q5)))))*SIGN(l.,X) 

415  ENDIF 

416  RETURN 

417  END 
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ADDITIONAL  PLOTS  OF  CURRENT  AND  RMS  ERROR 
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